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SUMMARY 

During  the  course  of  this  work,  we  have  developed  a  new  numerical 
method  for  the  integration  of  the  unaveraged,  time-dependent,  high-Reynolds 
number  Navier-Stokes  equations  governing  a  reacting  flow.  This  method, 
which  we  called  the  transport  element,  is  a  grid-free,  Lagrangian  field 
method  which  as  been  developed  for  the  simulation  of  reacting  flow,  and  to 
describe  mechanisms  of  shear  flow-combustion  interaction  which  have  been 
revealed  using  these  methods. 
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1.  Gagnon,  Yves,  Numerical  Investigation  of  Recirculating  Flow  at  Moderate 
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2.  Ng,  K.K.,  Vortex  Simulation  of  a  Confined  and  Perturbed  Mixing  Layer, 
M.Sc.  thesis,  February  1986. 

3.  Heidarinejad,  Ghassem,  Numerical  Simulation  of  a  Reacting  Shear  Layer 
Using  The  Transport  Element  Method,  Ph.D.  thesis,  February  1989? 

4.  Najm,  Habib,  Numerical  Investigation  of  the  Instability  of  Premixed 
Dump  Combustors,  Ph.D.  thesis,  June  1985. 
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2.  Loprenzo,  T. ,  Implementation  of  Reduced  Chemical  Kinetics  Schemes  in 
Vortex  Simulation  of  shear  Flow,  M.S  *-hesis,  started  September  1989. 

3.  Soteriou,  Marios,  Compressible  Vo  x  Method  for  the  simulation  of  high 
Mach  Number  Combustion,  Ph.D.  thesis,  started  September  1988. 
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2.  Krishnan,  A.  and  Ghoniem,  A.F.  "Origin  and  Manifestation  of  Flow- 
Combustion  Interaction  in  A  Premixed  Shear  Layer,"  Proceedings  of  the  22nd 
Symposium  (International)  on  Combustion,  the  Combustion  Institute, 
Pittsburgh,  PA,  pp.  665-675. (*) 

3.  Knio,  O.M.  and  Ghoniem,  A.F.  "Numerical  Study  of  A  Three-dimensional 
vortex  method,"  J .  Comput .  Phys . ,  in  press,  1989. 

4.  Heidarinejad,  G.  and  Ghoniem,  A.F.,  "Vortex  Simulation  of  the  Reacting 
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Aerospace  Sciences  Meetinq,  Reno,  Nevada,  January  9-12,  1989,  AIAA-89- 
0573.1*7 -  - 
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Workshop  on  The  Physics  of  Compressible  Turbulent  Mixing,  Princeton 
University,  Princeton,  N.J.,  October  24-57,  1988,  in  print. 

7.  Krishnan,  A.  and  Ghoniem,  A.F.,  "Numerical  Simulation  of  the  Structure 
of  A  Momentum/Gravity  Driven  Diffusion  Flame,"  27th  AIAA  Aerospace  Sciences 
Meeting,  Reno,  Nevada,  January  9-12,  1989,  AIAA-89-0485. ( * ) 

8.  Ghoniem,  A.F.,  Vortex  Methods  in  Turbulent  Reacting  Flow,  in  Numerical 
Approaches  to  Combustion  Modeling,  ed  by  E.  S.  Oran  and  J.  P.  Boris,  to  be 
published  by  the  AIAA,  1989. 

9.  Knio,  O.M.  and  Ghoniem,  A.F.  "Three-dimensional  Simulation  of  the 
Entrainment  Augmentation  Due  to  Streamwise  Vortex  Structures,"  27th  AIA& 
Aerospace  Sciences  Meeting,  Reno,  Nevada,  January  9-12,  1989,  AIAA-89-0574. 

10.  Knio,  O.M. ,  and  Ghoniem,  A.F.,  "Three  Dimensional  Vortex  Simulation  of 
Roll-up  and  Mixing  in  Shear  Flow,"  submitted  for  publication  at  the  Journal 
of  Computational  Physics,  June  1989. 


11.  Krishnan,  A.  and  Ghoniem,  A.F.,  "Simulation  of  the  rollup  and  mixing  in 
Rayleigh-Taylor  flow  using  the  vortex/transport  element  method,"  submitted 
for  publication  at  the  Journal  of  Computational  Physics,  August  1989. 

12.  Ghoniem,  A.F.  and  Krishnan,  A.  "Vorticity-combustion  interactions  in  a 
turbulent  reacting  jet,"  presented  at  The  12th  International  Colloquium  on 
Dynamics  of  Explosions  and  Reactive  Systems,  July  23-28,  1989,  Ann  Arbor, 
Michigan,  Proceedings  in  print. 

13.  Najm,  H.  and  Ghoniem,  A.F.,  "Flame-Vorticity-Acoustic  Interactions 
Leading  to  Combustor  Instability,"  the  AIAA/SAE/ASME/ASEE  25th  Joint 
Propulsion  Meeting,  Monterey,  CA,  July  10-12,  1989. (*) 
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14.  Ghoniem,  a.F.  and  Heidarinejad,  G.,  "The  structure  of  the  reaction  zone 


in  a  reacting  shear  layer;"  presented  at  The  12th  International  Colloquium 
on  Dynamics  of  Explosions  and  Reactive  Systems,  July  23-28,  1989,  Ann  Arbor, 
Michigan,  Proceedings  in  print. 


15.  Ghoniem,  A.F.,  "Numerical  Simulation  of  Turbulent  Combustion  using 
Vortex  Methods,"  Invited  five-lecture  series  at  the  University  of  Rome, 
Italy,  June  24-28,  1989. 


The  P.I.  delivered  a  number  of  invited  seminars  and  presentations.  A 
partial  list  is  presented  below: 


1.  Sandia  Meeting  on  Turbulent  Reacting  Flow,  General  Electric  Co., 
Schenectady,  October  1988. 

2.  DOD  and  EPA  Tyndall  Conference  on  Halon,  the  Ozone  Layer  and  Research 
on  Alternative  Chemicals,  Tyndall,  Florida,  November  1988. 

3.  university  of  Connecticut,  Storr,  November  1988. 

4.  university  of  California,  Berkeley,  April  1989. 

5.  Yale  University,  April  1989 

6.  Engineering  Foundation  Conference  on  Fluid  Mechanics  of  Engine 

Combustion,  Santa  Barbara,  May  1989 

7.  Gas  Research  Institute  Meeting  on  Fluid  Mechanics  of  Gas  Burners,  May 
1989. 

8.  Berkeley  Workshop  on  Vortex  Methods,  May  1989. 


(*)  Reprints  are  included  with  the  report. 
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INTERACTIONS  WITH  INDUSTRIAL  AND  GOVERNMENT  LABORATORIES  DURING  1988-1989: 


During  the  course  of  last  year,  we  have  started  and/or  continued 

collaborative  working  relations  with  the  following  industrial  or  governmental 

laboratories: 

1.  Wright-Patterson  Laboratory;  with  Dr.  M.  Roquemore  on  the  modeling  of  the 
jet  diffusion  flame. 

2.  General  Electric  Research  Center;  with  Dr.  Sanjay  Corea  on  the  study  of 
turbulent  premixed  flames  and  their  instability. 

3.  Sandia  National  Laboratory;  with  Drs.  R.  Lucht  and  John  Kelly  and  their 
associates  on  the  study  of  bluff-body  diffusion  flames. 

4.  Gas  Research  Institute;  with  Dr.  J.  Kezerle,  on  the  numerical  simulation  of 
bluff  body  flame  burners. 

5.  Ford  Motor  Company;  with  Mr.  C.  Kent  and  L.  Ramai,  on  numerical  simulation 
of  flame  propagation  in  enclosures  and  flow  during  intake  processes. 

6.  Shell  Oil  Co.;  with  Dr.  J.-C.  Ginestra,  on  numerical  simulation  of  flow, 
mixing  and  combustion  in  well-stirred  reactors. 

7.  Allison  Gas  Turbine;  with  Dr.  H.  Mongia,  on  numerical  simulation  of 
streamwibe  vortex  structures  in  thrust  augmentor  sections. 

8.  National  Institute  of  Standards  and  Technology;  with  Dr.  H.  Baum,  on 
numerical  simulation  of  uncontrolled  combustion. 

9.  Army  Atmospheric  Research  Laboratory;  with  Mr.  R.  Meyers,  numerical 
simulations  of  flow  over  complex  terrains. 
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TECHNICAL  REPORT 

VORTEX  SIMULATION  OF  REACTING  SHEAR  FLOW 

I.  BACKGROUND. 

Consider  the  class  of  unsteady,  high-Reynolds  number,  reacting-shear 
flows,  including  shear  layers,  jets,  wakes  and  recirculating  flows.  These 
flows  are  characterized  by  a  rapid  growth  of  successive  natural 
instabilities  which  lead  to  substantial  changes  in  the  geometry  of  the 
streamlines.  The  simulation  of  these  flows,  i.e.,  the  solution  of  the 
unsteady,  unaveraged  governing  equations,  requires  computational  schemes 
which  accommodate  these  geometrical  changes,  and  in  which  numerical 
diffusion  is  maintained  below  physical  diffusion  without  sacrificing 
numerical  stability.  One  way  to  achieve  this  is  to  utilize  schemes  based  on 
the  Lagrangian  formulation  of  the  conservation  equations. 

In  Lagrangian  methods,  grid  points,  where  the  flow  variables  are 
computed,  are  transported  along  particle  trajectories.  These  methods  offer 
natural  ways  to:  (1)  accommodate  the  severe  distortion  experienced  by  high- 
Reynolds  number  flows  as  they  evolve  through  different  states;  (2)  reduce 
the  artificial  diffusion  encountered  in  numerical  schemes  in  which  spatial 
derivatives  are  discretized;  and  (3)  account,  in  a  balanced  way,  for 
physical  diffusion  in  regions  where  the  two  modes  of  transport;  convection 
and  diffusion,  are  of  the  same  order  of  magnitude.  A  brief  discussion  of 
these  issues  is  presented. 

I.l.  NUMERICAL  ISSUES 

A  challenge  which  is  often  encountered  in  numerical  simulation  of 
reacting  shear  flow  is  how  to  capture  the  severe  distortion  of  the  flow  map 
which  results  from  the  nonlinear  growth  of  natural  flow  instabilities.  The 
saturation  of  these  instabilities  in  the  primary  flow,  which  can  normally  be 
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characterized  by  almost  parallel  streamlines,  results  in  a  secondary  flow 
with  strongly  curved  streamlines.  The  later  often  possesses  secondary 
instabilities  which  evolve  into  tertiary  flows.  To  accurately  capture  these 
changes  using  a  numerical  simulation,  a  very  large  number  of  fixed  grid 
points,  or  a  moving  grid  in  which  mesh  points  follow  the  distortion  of  the 
flow  field  is  needed.  The  latter  class  belongs  to  Lagrangian  schemes. 

Convection  is  the  dominant  mechanism  of  transport  at  high  Reynolds 
number.  Convection  in  the  cross-stream  direction,  perpendicular  to  the 
streamwise  direction,  is  increased  substantially  by  the  growth  of  natural 
flow  instabilities  and  their  saturation  into  fully-developed  flows. 
Meanwhile,  diffusive  transport  cannot  be  ignored.  Molecular  diffusion  is 
responsible  for  the  transport  of  vorticity  from  solid  walls  into  the 
interior  in  boundary  layers,  the  mixing  of  species  across  material  surfaces, 
and  the  transport  of  heat  and  species  across  laminar  flames.  Accurate 
simulation  of  diffusion  is  particularly  crucial  in  reacting  flows  where 
mixing,  most  often,  determines  the  burning  rate.  On  the  other  hand, 
excessive  numerical  diffusion  in  an  algorithm  can  artificially  stabilize  the 
flow.  Numerical  diffusion  acts  as  a  fictitious  source  of  molecular 
diffusivity  that  reduces  the  effective  Reynolds  number  of  the  computed  flow. 
In  combustion  calculations,  excessive  mixing  due  to  numerical  diffusion 
tends  to  increase  the  burning  rate.  To  capture  the  growth  of  flow 
instabilities  while  limiting  diffusion  to  practically  interesting  values  one 
needs  to  simulate  flows  at  Reynolds  number  0( 10^-10^). 

1.2.  CLASSIFICATION 

One  class  of  Lagrangian  methods,  which  has  successfully  been  used  in 
gas  dynamics,  utilizes  grids  to  discretize  flow  derivatives.,  e.g., 
Lagrangian  finite-difference  methods.  In  this  method,  grid  points  are 
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transported  along  particle  trajectories.  Shear  flow  can  lead  to  strong 
dictortion  of  the  grid  and  a  concomitant  loss  of  discretization  accuracy. 
To  maintain  accuracy  a  long  time  after  the  action  of  strong  shear,  mesh 
regularization  and  remeshing  become  important  [1-3]. 

On  the  other  hand,  grid-free,  Lagrangian  field  methods,  of  the  type 
described  in  this  report,  do  not  use  approximations  to  spatial  derivatives 
on  a  nonuniform  mesh.  Instead,  computational  elements  are  used  to  transport 
finite  values  of  the  spatial  gradients  of  the  variables,  specifically 
vorticity  and  scalar  gradients  (scalars  are  temperature  and  species 
concentrations).  Primitive  variables,  such  as  velocity  and  scalar 
concentrations,  are  obtained  by  integration  over  the  strength  of  the 
transport  elements.  These  method  are  labelled  field  methods  since  each 
element  induces  a  field  of  both  its  strength,  e.g.,  the  vorticity  of  an 
element  extends  over  a  finite  area,  and  of  the  primitive  variable  it  is 
transporting,  e.g.,  the  velocity  induced  by  a  vortex  element  extends  to  an 
infinitely  large  distance. 

The  goals  of  this  report  are  to  review  Lagrangian  field  methods  which 
have  been  developed  for  the  simulation  of  compressible  reacting  flow,  and  to 
describe  the  mechanisms  of  shear  flow-combustion  interaction  which  have  been 
revealed  using  these  methods.  These  goals  are  achieved  simultaneously  by 
progressively  introducing  more  complicated  models,  describing  the  necessary 
numerical  algorithms  and  presenting  their  results  in  a  form  most  relevant  to 
the  study  of  flow-combustion  interaction.  The  models  are  gradually  expanded 
until  they  include  many  of  the  physically  interesting  processes  in 
combustion. 


1.3.  ORGANIZATION 
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Vortex  methods,  a  particular  class  of  grid-free  Lagrangian  field 
methods,  have  been  used  to  obtain  solutions  of  the  momentum  equation.  These 
methods  are  based  on  the  discretization  of  the  vorticity  among  elements  of 
finite  area,  and  the  transport  of  these  elements  along  particle 
trajectories.  The  fact  that  vorticity  is  conserved  along  the  particle 
trajectory  in  a  two-dimensional,  uniform-density  flow  makes  these  methods 
particularly  simple.  However,  in  Section  II,  we  show  that  maintaining 
accuracy  requires  the  application  of  elaborate  vorticity-updating  schemes  as 
vortex  elements  are  moved  along  particle  trajectories  when  shear,  or  a 
strong  strain  field  is  present.  The  extension  and  application  of  vortex 
methods  to  three  dimensional  flows,  where  the  conservation  of  vorticity 
along  particle  trajectories  is  not  satisfied,  also  require  the  careful 
application  of  schemes  to  implement  the  effect  of  vortex  stretch  on  the 
strength  of  the  elements,  as  discussed  in  Section  III.  Solutions  using  the 
two  and  three-dimensional  methods  are  discussed  to  illustrate  some  of  the 
most  common  instabilities  encountered  in  nonreacting  and  reacting  shear 
flows  and  to  reveal  the  mechanisms  by  which  the  maturation  of  these 
instabilities  enhance  mixing,  and  hence  burning  in  a  reacting  flow. 

The  application  of  vortex  methods  to  nonuniform  density,  reacting  and 
compressible  flows  requires  compatible  Lagrangian,  grid-free  field  schemes 
to  compute  the  transport  of  scalars.  For  this  purpose,  we  developed  the 
transport  element  method  to  solve  the  convective-diffusive  scalar  transport 
equation.  Even  in  a  nonuniform  density,  incompressible  flow,  the  transport 
of  a  scalar  is  required  for  the  computations  of  the  density  field  which 
affects  the  flow  dynamically  via  the  generation  of  baroclinic  vorticity. 
The  formulation  of  this  method  is  described  in  Section  IV,  and  the  results 
of  its  application  to  compute  scalar  mixing  in  a  shear  layer  are  reviewed. 


10 


The  transport  element  method  is  then  combined  with  the  vortex  method  to 
solve  the  problem  of  nonuniform  density  shear  flow  in  Section  V.  Baroclinic 
vorticity  generation  is  one  of  the  important  mechanisms  by  which  combustion 
affects  the  fluid  dynamics  and,  thus,  the  computational  results  on  a 
nonuniform-density  shear  layer  are  reviewed  in  some  detail  in  this  section. 

Another  important  mechanism  of  flow-combustion  interaction,  namely 
reaction  extinction  due  to  the  formation  of  localized  regions  of  strong 
strains  as  instabilities  grow  into  their  nonlinear  range,  is  revealed  by  the 
results  of  incompressible  reacting  flow  models.  These  results  are  discussed 
in  section  VI,  after  reviewing  the  formulation  of  the  low-Mach  number 
combustion  model  and  the  extension  of  the  transport  element  method  to 
accommodate  the  reaction  terms  in  the  conservation  equations.  Three- 
dimensional  reacting  shear  layer  simulations  are  presented  to  illustrate  the 
complexity  of  the  mixing  pattern  encountered  in  these  flow  and  how  they 
affect  the  structure  of  the  burning  zone.  Compressible  reacting  flows 
exhibit  another  mechanism  of  combust ion- flow  interaction,  namely  the 
volumetric  expansion  associated  with  energy  release  within  the  reaction 
zone.  Results  of  computations  of  combustion  in  heterogeneous  and  homogenous 
compressible  shear  layers,  using  the  transport  element  method,  are  reviewed 
to  illustrate  the  origin  and  role  of  this  mechanism. 

The  bulk  of  this  report  is  devoted  to  a  brief  description  of  Lagrangian 
methods  for  the  simulation  of  reacting  flows  in  free  shear  layers,  and  a 
review  of  the  most  important  mechanisms  of  flow-combustion  interactions 
revealed  by  the  application  of  these  methods.  Section  VII  briefly  describes 
some  extension  of  vortex  methods  to  confined  shear  flow.  We  close  the 
report  with  some  thoughts  on  the  needs  for  future  development  in  the 
methodology  and  some  areas  of  application  which  have  not  been  explored. 
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II.  THE  VORTEX  ELEMENT  METHOD  IN  TWO  DIMENSIONS. 

The  vorticity  transport  equation  of  a  two-dimensu  .al,  incompressible, 
inviscid  flow  is: 

+  u  •  7w  =  0,  (1) 

where  Vxu  *  ■  and  V*u  -  0.  In  the  above,  u  *  (u,v)  is  the  velocity,  w  is 

the  vorticity,  x  »  (x,y),  t  is  time  and  V  -  (3/3x,3/3y).  If  x(X,t) 

describes  a  particle  path,  where  X  is  the  Lagrangian  coordinate  of  x  so  that 

X(X,0)  -  X,  then  Eq.  (1)  states  that  «<x(X,t),t)  -  w(X,0) ,i.e. ,  vorticity  is 

constant  along  particle  trajectory.  Thus,  if  the  vorticity  field  at  time  t 

-  0  is  divided  among  elements  which  move  along  particle  trajectories,  the 

strength,  i.e.,  the  vorticity  of  each  element  will  remain  the  same. 

2 

Moreover,  it  can  be  shown  that  u(x,t)  «  J  K(x-x')  w(x')  dx' ,  K(x)  -  -  l/2nr 
2  2  2 

(-y,x)  and  r  -  x  +y  ,  which  is  the  Biot-Savart  law.  This  Lagrangian 
formulation  of  vorticity  transport  is  the  basis  of  vortex  methods  [4,5]. 

II. 1.  NUMERICAL  SCHEME 

In  vortex  methods,  the  vorticity  field  is  discretized  into  a  number  of 
vortex  elements  of  finite  and  overlapping  cores: 

N  2 

w(x,t)  =  E  oo.  IT  f  i(x-X:  (X.  , t) ) ,  (2) 

i-1  1  6  ^  1 

vhere  is  the  vorticity  of  an  element,  N  is  the  total  number  of  vortex 

elements,  h  is  the  average  distance  between  the  centers  of  neighboring 

2 

elements  in  two  principal  directions,  h  *  h  h  ,  S  is  the  core  radius  of  a 

x  y 

2 

vortex  element,  and  f^  *  1/6  f(r/6)  is  the  core  function  describing  the 

distribution  of  vorticity  associated  with  an  element. 
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The  velocity  field  of  a  distribution  of  vortex  elements  is  obtained  by 
substituting  Eq.  (2)  into  the  Biot-Savart  law  and  performing  the 
integration.  The  result  is: 


u(x,t) 


•I  w.  h2  Ks(x-Xi(X,,t)), 


(3) 


where 

dXi  *  u()MX.  ,t),t),  (4) 

ar  ^  1 

while  Kj(x)  -  K(x)  K(r/S)  and  K(r)  -  2n  qJc  f(r')  r'  dr'. 

Equation  (2)  is  equivalent  to  expanding  a  function  w(x,t)  in  terms  of  a 

2 

number,  N,  of  kernel  functions,  f^,  located  at  X^  and  with  weights  w^h  . 
The  accuracy  of  the  discretization  depends  on  the  choice  of  f,  the  initial 

distribution  of  the  particles,  the  evaluation  of  the  values  of  i  -  1,  2, 

...  ,  N,  and  the  ratio  of  5/h  [6,7].  The  selection  of  a  core  function  to 
achieve  given  accuracy  was  extensively  discussed  in  a  number  of  theoretical 
analyses  [8-11].  For  the  evaluation  of  the  initial  values  of  we  found 
that  collocation  on  a  uniform  grid  provides  the  best  long  time  accuracy 
(collocation  on  a  nonuniform  grid  may  be  a  better  choice  when  the  initial 

vorticity  field  is  highly  nonuniform).  We  also  found,  using  extensive 

numerical  experimentation,  that  accurate  discretization  and  long-time 
accuracy  of  the  computed  flow  field  require  that  5  =  1.1-1. 5  h.  This  choice 
of  5/h  allows  for  strong  overlap  between  the  fields  of  the  vortex  elements. 
Thus,  the  local  value  of  vorticity  is  determined  by  the  contributions  of 
many  surrounding  elements. 

A  strong  strain,  associated  with  the  growth  of  perturbations  into  the 
nonlinear  stages  of  the  underlying  instability,  increases  the  distance 
between  neighboring  elements,  5x,  beyond  the  desired  maximum  value  of  h. 
Thus,  the  accuracy  of  the  spatial  discretization,  which  is  governed  by  5/h, 
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is  negatively  affected  and  unorganized,  random  motion  on  the  scale  of  h  is 
encountered.  This  deterioration  of  accuracy  is  connected  with  the  failure 
of  vortex  elements  of  fixed  cores  to  accurately  capture  the  distortion  of 
the  vorticity  field  locally.  To  circumvent  this  problem,  more  elements  are 
introduced  in  areas  where  5x  >  @h,  where  0  -  1.5.  The  circulation  of  the 
original  elements  is  locally  redistributed  among  the  newly  introduced 
elements,  and  the  local  value  of  vorticity  is  kept  constant. 

In  the  redistribution  algorithm,  conservation  of  vorticity  is  satisfied 
2  2 

by  dividing  h  and  5  such  that  5/h  remains  constant,  i.e.  the  core  radius 
of  an  element  is  decreased  as  the  element  is  exposed  to  strong  positive 
strain.  An  alternative  way  to  capture  the  effect  of  strain  on  the  geometry 
of  the  computational  elements  is  to  replace  a  circular  elements  by  am 
elliptical  element,  while  preserving  its  area,  with  its  major  axis  aligned 
with  the  direction  of  maximum  positive  strain.  If  this  was  done,  a  circular 
vortex  element  would  eventually  become  a  vortex  sheet  with  a  velocity  jump 
across  its  length  equal  to  the  local  value  of  the  vorticity  c*^.  For 
computational  convenience,  however,  we  replace  a  strained  circular  element 
by  several  circular  elements,  aligned  along  the  major  axis  of  strain,  but 
with  smaller  core  radii  than  the  original  element  [6]. 

II. 2.  PRIMARY  SHEAR  FLOW  INSTABILITIES 

We  used  shear  layers  as  test  cases  for  the  validation  of  the  numerical 
methods,  and  as  generic  problems  for  the  study  of  flow-combustion 
interactions.  We  start  with  a  two-dimensional,  incompressible,  nonreacting 
flow,  and  build  up  to  three-dimensional  flow,  variable  density  flow,  and 
compressible  reacting  flow. 

The  growth  of  a  sinusoidal  perturbation  on  an  infinite  shear  layer  is 
shown  in  figure  la  in  terms  of  all  the  vortex  elements  used  in  the 
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simulation  and  their  velocity  vector.  Figure  la  illustrates  the  adaptivity 
of  the  scheme;  more  elements  are  introduced  to  capture  the  straining  layers. 
The  growth  of  small  perturbations  in  a  spatially  developing  shear  layer  in 
which  the  flow  was  assumed  to  be  semi-infinite  is  depicted  figure  2  in  terms 
of  the  vortex  element  distribution.  In  both  the  temporal  and  spatial  shear 
layers,  numerical  results  for  the  growth  rate  of  small  perturbations  as  a 
function  of  the  perturbation  wavelength  were  found  to  agree  with  the  results 
of  the  linear  stability  theory.  The  nonlinear  regime  of  these  Kelvin- 
Helmholtz  instabilities  is  characterized  by  the  formation  of  large-scale 
structures  due  to  the  roll-up  of  the  vorticity  layer,  as  shown  in  figures  1 
and  2.  Results  in  this  regime  are  checked  for  convergence  and  against 
experimental  data.  Convergence  is  tested  by  establishing  the  grid 
independence  of  the  solution,  and  whether  the  solution  satisfies  certain 
differential  and  integral  constraints,  such  as  the  conservation  of  vorticity 
along  a  particle  path,  dw/dt  -  0. 

Figures  1  and  2  show  that  beyond  the  linear  range  of  the  instability, 
where  the  growth  is  exponential,  the  mean  growth  rate  of  the  shear  layer 
reaches  a  constant  value,  in  agreement  with  experimental  data.  Within  this 
range,  the  computed  results  were  used  to  evaluate  the  flow  statistics. 
Figure  3  shows  the  close  agreement  obtained  between  the  computed  results  and 
the  experimental  measurements.  We  note  that  mean  velocity  profiles  reach 
self-similarity  earlier  than  mean  fluctuations.  The  the  response  of  the 
layer  to  time-dependent  boundary  conditions  was  analyzed  by  modulating  the 
inlet  flow  at  frequencies  different  from  the  natural  shedding  frequency. 
Samples  of  the  results,  in  which  the  flow  is  forced  at  the  fundamental 
alone,  and  at  the  fundamental  and  the  subharmonic  at  the  same  time,  are 
shown  in  Figures  2b  and  2c,  respectively.  The  shear  layer  growth,  and  the 
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accompanying  rate  of  mixing,  is  enhanced  by  inducing  an  earlier  shedding, 
and  pairing  in  the  second  case,  through  the  application  of  external  forcing 
[7].  Multiple  pairing  can  be  achieved  by  lower  frequency  forcing,  as 
discussed  in  [13]. 

In  a  two-dimensional  flow,  the  source  of  fluctuation  is  the  formation, 
growth,  and  pairing  of  the  large-scale  vortex  structures  due  to  the  natural 
flow  instability.  Two  instabilities  are  observed  in  figure  2,  the  roll-up 
instability  which  leads  to  the  formation  of  large  eddies,  and  t_he  fairing 
instability  which  is  responsible  for  the  amalgamation  of  these  eddies 
downstream.  We  call  these  the  primary  instabilities  since  they  grow  in  the 
primary,  streamwise  direction.  Forcing,  which  can  be  used  to  either  promote 
or  suppress  these  instabilities,  was  shown  to  have  a  direct  impact  on  the 
values  and  signs  of  these  fluctuations,  suggesting  that  by  employing 
carefully  designed  forcing  functions,  one  can  control  the  interactions 
between  the  mean  flow  and  the  shear-layer  flow  [13].  Computations  of  the 
initial  stages  of  a  two-dimensional  turbulent  jet  were  presented  in  [14], 
III.  THE  VORTEX  ELEMENT  METHOD  IN  THREE  DIMENSIONS 

The  vorticity  transport  equation  in  an  incompressible,  three- 
dimensional,  inviscid  flow  is: 

~  +  u  •  Vu>  »  <i)  •  Vu  ( 5 ) 

where  oj  =  Vxu,  and  V  -  (3/3x, 3/3y, 3/3z) .  In  this  case,  u  *  (u,v,w)  and  x  * 
(x,y,z).  The  Lagrangian  form  of  Eq.  (5)  is  dw(X(X,t) ,t)/dt  = 
«( X(X,t)  ,t)  •  Vu,  where  Vu  is  the  strain  tensor  du^/dXy  An  equivalent 
expression  which  can  be  used  to  determine  the  vorticity  directly  is  the 
Helmholtz  theorem:  <o(X(X,t)  ,t)  *  VX(X,t) *«(X,0) ,  where  VX  is  the  Jacobian  of 
the  flow  map,  3x^/3Xj.  Moreover,  it  can  be  shown  that  u(x,t)  =  J  K(x-x')  x 
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c a(x')  dx',  where  K(x)  -  -1/An  x/r3,  which  is  the  Biot-Savart  law  in  three 

dimensions.  This  Lagrangian  formulation  is  the  basis  for  the  construction 
of  three-dimensional  vortex  methods  [15,16]. 

III.l.  NUMERICAL  SCHEME 

Vorticity  is  discretized  among  volume  elements,  of  side  h,  initially 
centered  around  X^,  by  collocation.  The  vortex  elements  are  then  moved 
along  particle  trajectories,  X^X^t),  while  their  vorticity  is  changed 
according  to  the  right-hand  side  of  Eq.  (5).  Thus, 

N  3 

«(x,t)  *  I  oj.  (t)  h*  f  c(x-X:  (X. , t) ) ,  (6) 

i=l  1  5  ^  1 

where  in  this  case,  f&(x)  -  1/S3  f(r/S)  and  the  rest  of  the  parameters  are 
defined  as  before.  Note  that,  here,  one  defines  strongly  overlapping  vortex 
balls  of  diameter  5,  and  that  the  core  function  is  spherically  symmetric, 
while  the  vorticity  vector  associated  with  an  element  is  u>^. 

The  total  vorticity  vector  of  an  element,  »j,h3,  is  expressed  more 

2 

naturally  in  terms  of  I\  and  51i  where  I\-  o^h  is  the  circulation  of  the 
element,  which  remains  constant  along  a  particle  path  in  an  inviscid  flow 
(Kelvin's  theorem),  and  61^  is  the  length  of  the  material  element  along  the 
vortex  line,  61^-  hi  t^/u^,  that  changes  as  the  material  lines  stretch 
(Helmholtz  theorem).  The  formulation  in  terms  of  (61,  T)  offers  a  natural 
regridding  method  which  is  used  when  61  >  0h  due  to  stretch.  In  this  case, 
a  vortex  element  is  divided  into  two  elements  along  the  vector  61,  while 
preserving  the  value  of  r.  The  conditions  necessary  for  the  accurate 
discretization  of  the  initial  vorticity  in  two  dimensions  are  valid  in  three 
dimensions,  e.g.,  5  >  h.  This  condition  must  be  satisfied  at  all  time. 
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The  velocity  field  induced  by  the  discrete  vorticity  distribution  is 
obtained  by  substituting  Eq.  (6)  into  the  Biot-Savart  law  and  integrating. 
The  results  can  be  written  as 

N 

U(x,t)  =  I  r.  81.  (t)  K.(x  -  X:(X.  ,t)),  (7) 

i=l  1  1  S  ^  1 

and 

dXi  -  u()C(X.  ,t),t),  (8) 

a r  ^  1 

51i(t)  -  \  <  W^+i'^  -  Vi(xi-i't))-  (9) 

where  K&(x)  -  K(x)  K(r/&)  and  K(r)  -  4n  QJr  f(r')  r'2  dr'.  To  reduce  the 
computations,  we  utilize  the  fact  that  vortex  lines  are  also  material  lines 
in  an  inviscid  flow.  Application  of  Eq.  (9)  requires  maintaining  data  on 
the  immediate  neighbors  in  the  direction  of  vorticity.  Thus,  one¬ 
dimensional  Lagrangian  grids  are  employed  to  describe  individual  vortex 
lines  as  arrays  of  vortex  elements  arranged  along  the  vortex  line.  The 
condition  that  a  vorticity  field  in  a  three-dimensional  free  space  be 
solenoidal,  V*<a  ■  0,  is  implicitly  satisfied  in  Eq.  (9).  Equations  (7)-(9) 
describe  the  vortex  filament  method  [17]. 

I I I. 2.  PROPAGATION  AND  INSTABILITY  OF  VORTEX  RINGS 

Before  describing  their  properties,  we  mention  that  vortex  rings,  in 
addition  to  providing  a  test  problem  for  the  numerical  scheme,  play  am 
important  role  in  combustion  theory  and  practice.  We  refer  the  reader  to 
the  experimental  literature  for  studies  on  flame-vortex-ring  interactions 
[18]  and  the  large-scale  structure  of  jet  diffusion  flames  [19]. 

It  has  been  shown  analytically  that  the  self-induced  velocity  of  a  thin 
vortex  ring,  a  «  R,  where  a  is  the  core  radius  and  R  is  the  ring  radius,  is 
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a  function  of  o/R  and  the  vorticity  distribution  within  the  core,  2(r/u) 

[20].  Numerically,  using  5  -  a,  the  dependence  of  the  self-induced  velocity 

on  <t/R  was  accurately  computed  when  overlap  was  maintained  between 

neighboring  vortex  elements  along  the  ring  axis  [21].  The  long-wave 

azimuthal  instability  of  a  thin  vortex  ring,  with  a  wavelength  X  >>  er,  was 

observed  when  the  ring  was  perturbed  along  its  axis  with  a  number  of  waves  n 

★  ★ 

-  2nR/X.  The  computed  unstable  wavenumber,  k  -  2n/n  ,  and  growth  rate  of 
these  waves  within  the  linear  range  agreed  with  the  prediction  of  the 
corresponding  linear  theory  [22].  The  growing  standing  waves  at  k*, 
contrary  to  the  spinning  stable  waves  at  all  other  wavenumbers,  expend  the 
flow  energy  in  stretching  the  waves  in  the  direction  opposite  to  the  ring 
self-induced  velocity. 

In  the  nonlinear  range,  the  growing  waves  form  almost  closed  loops  of 
vorticity  behind  the  original  ring,  as  shown  in  figure  4.  These  loops  are 
connected  to  the  original  ring  via  very  narrow  necks  that  can  be  pinched  off 
by  the  action  of  viscosity.  Each  loop  is  formed  of  two  vortex  rings  of 
opposite  signs  of  vorticity  separated  by  a  very  small  distance.  The 

separation  of  these  loops  from  the  parent  ring  may  lead  to  the  formation  of 
off-spring  vortex  rings  with  a  smaller  diameter  than  that  of  the  original 
ring.  This  would  result  in  an  interesting  cascading  to  smaller  scales;  a 
faster  decay  of  the  original  ring;  and  the  reduction  of  the  circulation  of 
the  parent  ring  [ 23 ] . 

The  study  was  extended  to  investigate  the  growth  of  short-wave 

instabilities,  X  -  a,  within  the  core  of  the  ring.  In  this  case,  one  must 

allow  the  core  of  the  ring  to  deform  under  the  action  of  the  growing 

perturbation.  Analysis  of  variations  within  the  core  requires  adequate 
resolution  of  its  vorticity  field  by  utilizing  vortex  elements  with  a  core 


19 


radius  8  <  <j,  i.e.,  several  elements  must  be  used  to  accurately  describe  the 

vorticity  field  within  the  core.  The  computed  value  k*,  shown  in  figure  is 

in  close  agreement  with  the  prediction  of  the  linear  theory  for  short- 

wavelength  instability  in  deformable  vortex  rings.  The  value  of  k*  depends 

on  the  vorticity  distribution  within  the  ring  core,  2(r/o).  figure  5  shows 

★ 

that  the  value  of  k  ,  predicted  from  the  short-wave  analysis,  is  closer  to 
the  experimental  data  than  that  predicted  by  the  previous  long-wavelength. 

Spectral  analysis  of  the  field  of  an  unstable  ring  shows  that  as  the 
fundamental  instability  reaches  saturation,  its  harmonic  becomes  unstable 
and  starts  to  grow.  The  mechanism  of  excitation  of  this  frequency, 
displayed  in  figure  6,  is  associated  with  the  formation  of  hairpin  vortices 
at  half  the  wavelength  of  the  original  perturbation.  Examining  the 
vorticity  field  after  the  saturation  of  the  first  instability  reveals  the 
presence  of  a  strong  streamwise  vorticity  component  with  alternating  signs 
as  one  moves  along  the  axis  of  the  original  ring  [24].  This  component  is 
the  result  of  the  tilting  of  the  original  vortex  lines  into  the  streamwise 
direction.  The  tilting  occurs  as  the  azimuthal  instability  grows.  Within 
each  azimuthal  wave,  two  vortices  of  opposite  signs  are  formed. 

The  scheme  was  also  used  to  study  the  three-dimensional  instability  of 
the  initial  stages  of  an  axisymmetric  jet  [17].  It  was  found  that  the  jet 
rolls  up  into  vortex  rings  which  experience  a  similar  instability  to  that  of 
an  isolated  vortex  ring. 

III. 3.  SECONDARY  INSTABILITIES  OF  SHEAR  FLOW. 

Beyond  the  primary  instability,  shear  flows  develop  secondary 
instabilities  which  leads  to  the  generation  of  streamwise  vorticity  (the 
vorticity  vector  points  in  the  streamwise  direction),  an  important  mechanism 
of  mixing  enhancement  [24-26].  To  simulate  this  process,  the  shear  layer  is 


20 


initially  perturbed  in  the  streamwise  and  spanwise  directions,  and  the 
vortex  scheme  is  modified  to  accommodate  the  strong  strain  field  which 
develops  in  the  plane  normal  to  the  initial  vorticity.  This  is  accomplished 
by  redistributing  the  vorticity  field  among  a  larger  number  of  elements 
arranged  in  the  direction  of  maximum  strain. 

In  the  transport  element  scheme,  each  vortex  element  is  defined  by  its 
circulation  and  two  material  vectors:  one  in  the  direction  of  vorticity  and 
the  other  in  the  direction  of  maximum  strain.  This  is  equivalent  to  using 
two-dimensional  Lagrangian  grids  to  describe  planes  of  constant  circulation. 
The  stretch  of  the  sides  of  the  grid  can  be  used  to  update  the  vorticity 
associated  with  each  element  in  the  plane.  Alternatively,  a  grid-free 
stretching  scheme  cam  be  constructed  on  the  basis  of  Eq.  (5).  In  this 
scheme,  the  term  a».Vu  is  computed  for  each  element  by  analytically 
differentiating  the  velocity  expression  in  Eq.  (7).  The  computational 
effort  is  almost  the  same  in  both  schemes,  however  the  bookkeeping  in  the 
first  scheme  is  greater  [26]. 

The  evolution  of  a  temporal  shear  layer  is  displayed  in  the  next  three 
figures.  Figure  7  depicts  the  distortion  of  the  material  surface  initially 
aligned  along  the  midsection  of  the  layer,  where  u  -  0,  in  terms  of  the  two- 
dimensional  grid  used  to  discretize  its  vorticity  (this  is  one  of  several 
such  surfaces  used  to  discretize  the  vorticity).  Figure  8  shows  the 
distribution  of  the  streamwise  vorticity  contours  on  a  y-z  streamwise  plane 
(a  streamwise  plane  is  normal  to  the  streamwise  direction)  which  cuts 
through  the  core  of  the  spanwise  structure  at  the  middle  of  the  domain. 
Figure  9  depicts  the  streamwise  vorticity  on  a  y-z  plane  which  cuts  through 
the  braids  of  the  spanwise  structures  (the  principal  axis  of  this  structure 
coincides  with  the  spanwise  direction).  The  distribution  of  the  streamwise 
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vorticity  indicates  that  during  the  growth  of  the  strearawise  perturbation, 
the  growth  of  the  spanwise  perturbation  is  suppressed,  and  the  flow  remains 
almost  two-dimensional.  The  growth  of  the  streamwise  instability  (not 
shown)  matches  that  of  a  two-dimensional  flow  until  the  time  the  streamwise 
instability  saturates,  t  -  8.  The  saturation  is  accompanied  by  the 
formation  of  a  spanwise,  c  ir  -;ally  shaped  vortex  core  whose  axis  is 
perpendicular  to  the  mean  flow  direction.  The  formation  of  this  core  is 
accompanied  by  a  secondary  flow  whose  streamlines  are  almost  circular. 

Beyond  the  saturation  of  the  streamwise  instability,  t  >  8,  the  growth 
rate  of  the  spanwise  instability  increases  and  the  vortex  core  starts  to 
deform.  This  deformation,  and  the  concomitant  tilting  of  the  vorticity 
vectors  from  the  spanwise  direction  into  the  streamwise  direction,  leads  to 
the  establishment  of  a  streamwise  vorticity  component  with  an  alternating 
sign  along  the  axis  of  the  spanwise  cylindrical  structures.  The  deformation 
of  the  spanwise  cylindrical  cores  is  known  as  the  translative  instability. 
The  streamwise  vorticity  associated  with  this  instability  is  shown  in  figure 
8.  The  fact  that  the  cores  are  covered  from  both  sides  by  two  rows  of 
streamwise  vortex  rods  is  explained  next. 

The  other  source  of  streamwise  vorticity  is  the  deformation  of  the 
vortex  lines  of  the  braids  in  the  cross-stream  direction,  as  shown  in  figure 
7.  The  vorticity  component  in  the  cross-stream  direction  is  transformed 
into  a  streamwise  component  by  the  action  of  the  velocity  gradient,  3u/3z, 
where  u  is  the  velocity  component  in  the  x-direction.  Along  the  spanwise 
direction,  streamwise  vorticity  changes  its  sign  every  half  wavelength. 
This  configuration  is  unstable,  and  each  half  wavelength  rolls  up  to  form  a 
vortex  rod  aligned  with  the  local  direction  of  the  braids,  as  shown  in 
figure  9.  The  vorticity  within  these  rods  is  constantly  amplified  as  the 
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strain  field  in  the  streamwise  direction,  generated  by  the  large  spanwise 
cores,  strains  the  flow  along  the  direction  of  the  braids.  With  elongation, 
these  rods  wrap  around  the  spanwise  cores  and  produce  the  distribution  shown 
in  figure  8. 

The  deformation  of  the  flow  field  due  to  the  evolution  of  the  two  and 
three-dimensional  instabilities  leads  to  a  substantial  enhancement  of  mixing 
within  the  large  structures  by  generating  strong  entrainment  currents  of 
fluids  from  both  sides  of  the  shear  layer  towards  the  vortex  center.  The 
mixing  enhancement,  and  its  effect  on  the  burning  rate  are  discussed  in 
detail  in  the  next  sections. 

IV.  THE  TRANSPORT  ELEMENT  METHOD 

Given  that  s  is  a  passive,  nondiffusive  scalar,  the  conservation 
equations  for  s,  and  g  -  Vs,  in  an  incompressible  flow  are: 


ds 

at = 


o. 


and 


g  x  <o. 


(10) 


(ID 


Thus,  s  remains  constant  along  a  particle  path,  while  g  changes  due  to  the 
straining  and  rotation  of  the  material  lines  by  the  local  strain  field  and 
vorticity.  If  the  material  is  exposed  to  a  strong  strain  in  the  direction 
normal  to  the  gradient,  the  value  of  g  must  increase  by  the  same  amount  as 
the  stretch  in  the  material  element.  This  can  be  seen  by  deriving  an 
explicit  equation  which  relates  the  changes  in  g  -  |g|  to  the  variation  of 
material  elements,  or  the  distortion  of  the  flow  map.  This  is  done  by 
expanding  Eq.  (11)  in  terms  of  g  n,  and  implementing  kinematical  relations 
that  describe  the  variations  of  n  =  g/g,  where  n  is  the  unit  vector  normal 
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to  an  isoscalar  line,  i.e.,  a  line  along  which  s  is  constant.  After  some 
lengthy  manipulations,  we  get: 

^3n  *  -g  (n*7u  +  j  n  x  u  -( 1* (n*7us) )1) ,  (12) 

where  Vus  is  the  symmetric  part  of  the  strain  tensor  Vu  and  1  is  the  unit 
vector  normal  to  n.  Moreover,  g  ■  (ds/dn)  n  -  (6s/5n)  n,  where  6s  is  the 
variation  of  s  across  a  small  material  line  6n.  In  two  dimensions  and  for 
an  incompressible  flow,  the  variation  of  a  material  vector  element  61  can  be 
shown  to  be  governed  by  a  similar  equation  to  Eq.  (12). 

IV. 1.  NUMERICAL  SCHEME. 

From  the  above  discussion  on  the  relationship  between  scalar  gradients 
and  the  deformation  of  the  flow  map,  it  follows  that  g/61  -  constant  along  a 
particle  path,  and  that  the  scalar  gradient  can  be  computed  from  the 
following  relations  [7]: 


N  2 

g(x,t)  -  L  g.(t)  h  f6(x  -  X^X^t)), 

where 

6s.  61. (t) 

9i(t)  *  72"  ni(t)' 


(13) 


(14) 


and  ^ ( Xi , t )  is,  as  before,  a  particle  path.  Equation  (13)  is  based  on  the 
expansion  of  g  in  terms  of  the  core  function  f^,  similar  to  Eq.  (2).  Since 
an  isoscalar  line  is  a  material  line  in  a  nondiffusive  field,  61^  can  be 
updated  as:  61^(t)  -  ( Xi+i~Xi_i )/2,  while  n^l^  *  0.  Thus,  it  suffices  to 
move  the  centers  of  the  transport  elements  while  remembering  the  near 
neighbors  at  t  -  0.  As  in  the  vortex  method,  when  the  distance  between 
neighboring  elements  in  the  direction  of  maximum  strain  exceeds  a  certain 
maximum,  an  element  is  inserted  between  two  neighboring  elements.  The  total 
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of  Sl^'s  for  the  two  original  elements  are  distributed  among  the  three 
2  2 

elements,  while  h  and  5  are  adjusted  so  that  the  total  material  area  is 
conserved,  keeping  Ss^  the  same. 

For  a  variable-density  flow,  the  above  analysis  is  modified  to  reflect 
the  fact  that,  in  this  case,  sin  equation  similar  to  Eq.  (12)  can  be  derived 
with  g  replaced  by  p  51,  and  the  expression  of  g^  in  Eq.  (13)  changes  to 
[12,27]: 


9i(t) 


6si  Sl.(t)  p.(t) 

pTTo)  h2 


n^t) . 


(15) 


The  value  of  p  is  computed  using  the  relation  p  T  -  constant  in  the  low  Mach 
number  approximation  (see  Section  V).  Given  the  location  and  strength  of 
the  transport  elements,  the  scalar  concentration  are  computed  by  direct 
integration  over  the  fields  of  the  transport  elements 

N  2 

s ( x, t )  «  .E  gi(t)  hi  •  TOs(x-xi(Xi,t)),  (16) 

2 

where,  in  two  dimensions,  VG^(x)  -  (x,y)/2nr  <(r/6),  and  in  three 

dimensions,  VG^(x)  -  (x,y,z)/4nr^  <(r/5).  This  formulation  is  fully 

compatible  with  the  vortex  method  since  all  the  information  needed  to 
compute  the  scalar  transport  are  already  part  of  the  vortex  computations, 
including  all  the  expressions  of  the  Green  functions.  For  extended 
derivations,  see  [24,27]. 

The  effect  of  molecular  diffusion  in  the  vortex  and  transport  element 

methods  can  be  modeled  by  expanding  the  cores  of  the  elements  according  to 

2  2 

the  following  relation,  5  (t+At)  =  5  (t)  +  2  a  At,  where  At  is  the  time  step 
and  a  is  the  molecular  diffusivity  [6].  This  relation  is  obtained  by  direct 
suostitution  of  Eq.  (11)  into  the  diffusion  equation.  A  limit  should  be 
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imposed  on  the  maximum  allowable  value  of  8  to  maintain  the  spatial  accuracy 
of  the  calculations.  Beyond  6.  an  element  should  be  subdivided  into  a 
number  of  smaller  elements.  Another  scheme  for  implementing  the  effect  of 
diffusion  without  expanding  the  cores  was  proposed  in  [28].  This  scheme  is 
used  in  conjunction  with  the  three-dimensional  transport  element  method  in 
Section  IV. 3. 

IV. 2.  SCALAR  MIXING  IN  SHEAR  LAYER 

The  transport  element  method  was  used  to  study  the  transport  of  species 

in  a  two-dimensional,  heterogeneous  shear  layer  [12].  Figure  lb  shows  the 

isoscalar  lines  of  the  shear  layer  of  figure  la  when  the  diffusivity  is 

zero.  Numerically,  we  found  that,  provided  the  field  is  accurately 

discretized  initially,  the  condition  ds/dt  =  0  is  satisfied  if  the  core 

radii  of  the  elements  are  allowed  to  decrease  at  the  rate  described  in 

2  2 

Section  II,  i.e.  if  h/8  remains  constant  as  the  elements  are  deformed. 
This  also  guarantees  that  the  scheme  can  capture,  without  introducing 
numerical  diffusion,  the  large  scalar  gradients  that  arise  from  the  strong 
deformations  in  the  flow  which  accompany  the  evolution  of  the  instability. 

Figure  10  shows  a  comparison  between  the  computed  and  the  measured 
values  of  the  mean  concentration  and  the  root-mean  square  of  the 
concentration  fluctuations  in  the  two-dimensional  spatially  developing 
mixing  layer  of  figure  2a  [12,29].  The  computational  results  were  obtained 
for  a  range  of  Peclet  number  between  10 ^  and  10^,  so  that  the  dominant 
transport  process  was  convection  and  not  diffusion.  However,  the  effect  of 
species  diffusion  was  incorporated  to  study  mixing.  The  effect  of  diffusion 
on  the  mean  scalar  distribution  is  very  small  since  the  overall 
concentration  field  is  established  by  the  convective  field,  also  called  the 
entrainment  currents.  Due  to  the  roll-up  of  the  vorticity  layer,  fluid  from 
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both  sides  of  the  shear  layer  are  engulfed  into  the  large  structure  and,  on 
the  average,  mean  values,  between  the  two  extremes,  can  be  encountered. 

The  root-mean  square  of  the  concentration  fluctuations  exhibit  strong 
dependence  on  the  Peclet  number.  Its  maximum  value,  0.5,  can  only  be 
achieved  at  very  high  Peclet  number  where  strong  unmixedness  is  present 
inside  the  cores  of  the  eddies.  As  the  effect  of  molecular  diffusion 
increases,  it  homogenizes  the  cores  and  the  fluctuations  drop  below  0.5.  In 
this  case,  the  profiles  show  a  mixed  region  of  constant  value  of 
fluctuation.  For  all  values  of  the  Peclet  number,  the  fluctuations  never 
reached  zero  inside  the  cores  indicating  that  the  fluid  does  not  reach 
complete  homogeneity.  Another  interesting  feature  of  these  mixing  flows,  is 
the  presence  of  more  high  speed  fluid  than  low  speed  fluid  inside  the  cores 
[12],  This  mixing  asymmetry,  which  is  a  direct  consequence  of  the  unequal 
velocities  across  the  interface  of  the  layer,  will  be  shown  to  play  an 
important  role  in  the  distribution  of  products. 

IV. 3.  ENTRAINMENT  ENHANCEMENT  DUE  TO  3D  INSTABILITIES 

The  three-dimensional  transport  element  method  is  based  on  the 
relationship  between  the  distortion  of  the  flow  map  and  the  scalar 
gradients,  as  in  the  two-dimensional  method.  In  3D,  the  distortion  of  the 
flow  is  represented  by  the  changes  in  the  magnitude  and  direction  of  a 
material  area  element  initially  aligned  with  an  isoscalar  surface  3A.  This 
relationship  was  derived  in  [24,30]  for  a  variable  density  flow; 

dtln  g  "  Htln  (p  (17) 

Equation  17  indicates  that  g/( p  6A)  =  constant  along  a  particle  path.  To 
implement  this  relationship  in  the  calculations  of  the  scalar  gradients,  one 
must  keep  track  of  the  area  of  the  elements  associated  with  computational 
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points  moving  along  particle  trajectories.  Similar  schemes  are  used  in  the 
vortex  element  method  to  monitor  the  effect  of  the  strain  field  on  the 
evolution  of  the  vorticity  field.  Thus,  the  computation  of  the  scalar 
transport  is  a  natural  extension  of  the  vortex  element  method  since  it  does 
not  require  additional  computational  effort. 

The  entrainment  due  to  the  development  of  three-dimensional 
instabilities  in  shear  flow  is  depicted  in  figures  8  and  9  where  we  show  the 
scalar  contours  on  the  same  planes  where  the  streamwise  vorticity  contours 
are  displayed.  The  extra  entrainment,  over  what  is  observed  in  two- 
dimensional  calculations,  is  induced  by  the  action  of  the  spanwise 
deformation  of  the  vortex  core  and  by  development  of  the  streamwise  vortex 
rods  within  the  braids,  as  shown  in  figures  8  and  9,  respectively.  Hie 
total  entrained  fluid  within  the  large  eddies  is  measured  by  the  size  of  the 
eddy  in  the  cross-stream  direction,  and  shown  in  figure  11  for  both  the  two 
and  the  three-dimensional  calculations  [24,26].  The  effect  of  this  extra 
entrainment  on  the  rate  of  burning  will  be  discussed  in  Section  VI. 2. 

V.  VORTEX  METHODS  FOR  VARIABLE-DENSITY  FLOWS 

In  combustion  phenomena  of  practical  interest,  the  flow  is  density 
stratified,  i.e.  finite  and  large  density  gradients  exist  between  material 
layers,  due  to  heat  release  and  the  variation  in  molecular  weight  of 
reacting  species.  The  effect  of  density  gradients  is  to  introduce  a 
baroclinic  source  term  in  the  vorticity  transport  equation.  In  this 
section,  we  show  how  this  term  is  included  in  the  simulation,  and  describe 
its  effect  on  shear  layers.  Considering  an  incompressible,  inviscid, 
variable  density,  the  conservation  equations  are: 


(18) 
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dw  1  _  _  _ 

^  =  -^7px7p  +  M.  9u, 

P 


(19) 


where  7p  »  -  p  du/dt  and  Vp  is  computed  using  the  transport  element  method 
as  described  before.  Extensions  to  gravity-driven  at  finite  Froude  number 
are  considered  in  (27,31].  To  take  the  extra  term  in  Eq.  (19)  into  account, 
the  computational  algorithm  proceeds  in  fractional  steps:  (1)  the  vorticity 
is  transported  without  change  along  the  particle  trajectory,  Eq.  (8);  (2) 
the  material  acceleration  along  the  particle  path  is  computed  to  find  the 
pressure  gradient;  (3)  the  change  in  the  density  gradient  along  the  particle 
path  is  evaluated  using  the  transport  element  method,  Eq.  (14);  and  (4)  the 
vorticity  is  updated  according  to  the  discrete  analog  of  Eq.  (19), 


dr. 

i 

ar 


Vp.  ,du> 

— * x  Wi 


(20) 


V.l.  BAROCLINIC  EFFECTS  IN  SHEAR  LAYER. 

Computations  of  a  variable  density  shear  layer  were  performed  to 
investigate  the  effect  of  baroclinic  vorticity  generation  on  the  growth  rate 
of  the  Kelvin-Helmholtz  instability,  the  rate  of  entrainment  and  the  rate  of 
expansion  of  a  shear  layer.  Figure  12  shows  a  sample  of  a  two-dimensional 
simulation  in  which  the  fast  top  stream  is  five  times  heavier  than  the  slow 
bottom  stream.  This  configuration  resembles  the  experimental  set-up  used  to 
study  a  shear  layer  between  premixed  reactants  moving  in  the  top  stream  and 
products  of  combustion  moving  in  the  bottom  stream  [ 34  ] .  The  computed 
initial  growth  rate  and  phase  velocity  of  the  growing  waves  compared  well 
with  the  results  of  the  linear  stability  theory  [12,  35). 

In  the  nonlinear  range,  the  computational  results  indicate  that  density 
variation:  (1)  induces  a  net  convective  motion  on  the  eddy  in  the  direction 
of  the  heavy  stream;  (2)  enforces  entrainment  asymmetry  on  the  growing  eddy 
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which  results  in  the  presence  of  more  light  fluid  than  heavy  fluid,  by 

volume,  within  the  core;  and  (3)  destabilizes  the  flow  leading  to  the 

formation  of  large  scale  structures  within  the  large  scale  structure.  It 

was  also  observed  that  baroclinic  vorticity  enhances  entrainment  over  that 

of  the  uniform-density  case.  Using  a  Galilean  transformation  between  the 

temporal  and  spatial  shear  layer,  i.e.,  dS/dx  -  1/U  dS/dt  where  u  is  the 

c  c 

finite  convective  velocity,  one  can  show  that  a  shear  layer  in  which  the 
heavy  stream  moves  faster  than  the  light  stream  will  grow  more  slowly  in  the 
streaunwise  direction  than  a  uniform  density  layer,  in  agreement  with 
experimental  observations  ( 36 ] . 

Similar  changes  were  found  in  the  results  of  a  three-dimensional 
simulation  of  a  variable  density  flow.  Figure  13  shows  a  sample  of  these 
results  for  the  case  of  a  top,  fast  stream  is  twice  as  heavy  as  the  bottom, 
slow  stream  [30].  The  results  are  represented  by  the  grid  used  to  describe 
the  material  surface  initially  aligned  with  the  midsection  of  the  shear 
layer,  as  in  figure  7.  The  motion  of  the  developing  eddy,  its  asymmetry 
with  respect  to  the  top  and  bottom  streams  and  some  added  asymmetry  in  the 
spanwise  direction  are  clearly  seen  in  the  figure.  Secondary  instabilities, 
leading  to  the  deformation  of  the  spanwise  core  and  the  formation  of 
streaunwise  vortex  rods  within  the  braids,  are  similar  to  those  encountered 
in  the  two-dimensional  simulations. 

The  dynamic  origin  of  these  phenomena  is  the  generation  of  vorticity 
due  to  the  baroclinic  torque.  The  misalignment  of  the  density  and  pressure 
gradients  within  the  growing  eddy  leads  to  the  generation  of  additional 
vorticity  with  two  opposite  signs  on  the  two  sides  of  the  shear  layer,  as 
shown  in  figure  12.  In  the  case  considered  here,  vorticity  of  a  similar 
sign  as  that  of  the  shear  layer  is  generated  on  the  light,  slow  fluid  side 
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while  vorticity  of  the  opposite  sign  is  generated  on  the  fast,  heavy  fluid 
side.  This  redistribution  of  vorticity  within  the  growing  eddy  leads  to 
more  entrainment  from  the  lower,  light  fluid  side  of  the  layer  and  a  net 
streamwise  motion  of  the  large  eddy  in  the  direction  of  the  heavy  fluid. 
The  computed  entrainment  ratio  and  finite  eddy  velocity  were  found  to  agree 
with  experimental  results  (32).  In  Section  VI. 4,  we  show  how  baroclinic 
vorticity  generation  affects  combustion  in  premixed  shear  layer. 

The  effect  of  baroclinic  vorticity  on  the  instability  of  a  jet  was 
presented  in  [14],  without  and  with  the  effect  of  Froude  number.  The  role 
of  baroclini'-  vorticity  in  determining  the  structure  of  low-speed  jet 
diffusion  flames  is  analyzed  in  [27].  The  extension  of  the  calculation  of  a 
variable  density  flow  to  a  premixed  shear  layer,  and  the  corresponding  role 
of  baroclinic  vorticity  are  discussed  in  Section  VI. 4. 

VI.  VORTEX  METHODS  FOR  REACTING  FLOWS 

The  transport  element  method  was  extended  to  reacting  flow  [7,35,37]. 
In  combustion  systems,  heat  release  at  low  Mach  number  leads  to  the 
generation  of  an  irrotational  velocity  field,  ?♦,  which  when  superimposed  on 
the  existing  vorticity-induced  rotational  velocity,  u^,  represents  the  total 
velocity  in  a  reacting  flow.  The  low  Mach  number  approximation  is  employed 
to  filter  out  the  pressure  waves  and  to  render  the  pressure  independent  of 
the  density  in  unconfined  flows  [38].  Employing  the  velocity  decomposition, 
the  governing  equations  of  combustion  can  be  written  as  follows: 

u  -  uw+  V$,  (21) 


72*  -  I  *  , 
*  T  Ht  ' 


=  -2  Vp  x  7p  -  a>  (V.u)  +  (w.V)  u  +  |j-  720), 
p  e 


(22) 


(23) 
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a?  *  r  v2t  +  Q  w  ' 

(24) 

e 

^ck  1  2 

ar*pL  7  ck  +  wk  ' 

e  e 

(25) 

p  T  *  constant, 

(26) 

where  T  is  temperature,  c  is  species  concentration,  W  is  the  rate  of 
creation/destruction  of  species,  Q  is  the  enthalpy  of  reaction,  Rg,  Pg  and 
Lg  are  the  Reynolds,  Peclet  and  Lewis  numbers,  respectively.  Here,  we 
assume  that  (1)  the  transport  properties  are  constant  except  for  mass 
diffusivity  which  is  replaced  by  a  density-averaged  value;  and  (2)  molecular 
weights  of  reactants  and  products  are  the  same.  The  pressure  gradient  is 
obtained  from  the  momentum  equation:  Vp/p  -  -  a  +  l/(pRg)  V  u,  where  a  ■ 
dxydt.  Both  velocity  components,  and  V<$>  must  satisfy  the  velocity 
boundary  condition  in  the  direction  normal  to  the  boundaries  of  the  domain. 
Using  the  principle  of’  vortex  decomposition,  i.e.,  discretizing  the  source 
term  in  Eq.  (22)  in  a  form  similar  to  Eq.  (2),  Green's  function  solution  of 
Eq.  (22)  can  be  written  as  17]: 


V*(x,t) 


N 

l 

i*l 


1 

T. 


(g£)i  h?( t)  7G5(x-Xi(X.,t)). 


(27) 


The  energy  equation  and  the  species  transport  equations,  Eq.  (24)  and 
(25),  respectively,  are  solved  using  the  transport  element  method  in  three 
fractional  steps:  convection,  diffusion  and  reaction.  The  reaction 
fractional  step  is  implemented  by  changing  the  strength  of  the  transport 
elements  according  to  the  following  expression  [7]: 
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where  k  is  the  total  number  of  species  and  sk+1  =  T. 

In  the  following,  we  review  four  solutions  of  reacting  shear  layers, 
computed  using  the  transport  element  methods.  The  first  three  solution  are 
obtained  for  the  reaction  between  two  streams  of  fuel  and  oxidizer,  and  the 
fourth  is  for  a  reaction  between  a  stream  of  premixed  reactants  and  a  stream 
of  products.  In  the  first  two  cases,  an  incompressible  flow  model  was  used 
to  eliminate  the  complexity  of  the  effect  of  heat  release  on  the  flow 
dynamics.  In  the  third  and  fourth  cases,  a  compressible  flow  model  was 
used.  In  all  cases,  we  focus  on  the  mechanisms  of  shear  flow-combustion 
interactions  under  different  physical  conditions. 

VI. 1.  THE  REACTING  SHEAR  LAYER 

A  computation  of  a  reacting  shear  layer  between  a  fuel  stream,  F,  and 
an  oxidizer  stream,  O,  in  a  spatially  developing  flow  was  described  in 
[39,40].  The  incompressible  flow  model,  p  «  constant,  was  used  to 
investigate  how  the  flow  which  evolves  due  to  the  roll-up  instability 
affects  the  combustion  process.  The  chemistry  is  modeled  by  one-step 
reaction  kinetics,  F  +  O  ->  P,  and  the  reaction  rate  is  taken  as  W  -  A^  cR 
c_  exp(-T/T  ),  where  is  the  frequency  factor  and  T  is  the  activation 
energy  normalized  with  respect  to  the  gas  constant.  A  sample  of  the 
results,  obtained  for  Ta/TR-  10,  A^-  8000  and  Q  »  1,  is  shown  in  figure  14, 
revealing  a  strong  similarity  between  the  temperature,  product  concentration 
and  the  vorticity  contours,  being  highest  within  the  cores  and  reaching  very 
small  values  within  the  braids.  A  large  eddy,  defined  by  a  closed  set  of 
vorticity  contours,  figure  14c,  can  also  be  described  by  a  set  of  closed 
product  concentration  contours,  around  the  same  origin  and  of  the  same  shape 
as  the  vorticity  contours,  figure  14b.  This  similarity  confirms  the 
important  role  of  advection,  which  is  governed  by  vorticity  in  this  case,  in 


determining  the  local  concentration  field,  the  rate  of  mixing  and  the  rate 
of  burning. 

The  effects  of  the  shear  layer  roll-up,  the  Damkohler  number,  Reynolds 
number  and  reactants  ratio  across  the  layer  are  shown  in  figure  15.  These 
results  were  obtained  for  a  temperature-independent  reaction  in  which  W  -  Dq 
cF  Cp,  where  Dg  is  the  Damkohler  number.  The  roll-up  of  the  shear  layer 
enhances  product  formation  over  the  laminar  layer  by  inducing  strong 
entrainment  currents  into  the  cores.  A  drop  in  the  total  product  formation 
in  the  early  stages  is  due  to  the  thinning  of  the  reaction  zone  by  the 
strain  field.  Product  formation  depends  strongly  on  the  Damkohler  number. 
This  dependence,  however,  vanishes  and  the  total  amount  of  product  formed 
within  the  layer  reaches  an  asymptotic  value  at  Damkohler  number  on  the 
order  of  20,  based  on  the  length  scale  of  the  large  structures.  This  is  the 
fast  chemistry  limit  at  which  the  rate  of  product  formation  depends  only  on 
the  mixing  rate. 

At  high  Reynolds  numbers,  of  an  order  of  10000  based  on  the  length 
scale  of  the  large  eddies,  product  formation  is  a  weak  function  of  the 
Reynolds  number.  This  indicates  that  while  mixing  is  strongly  enhanced  by 
convective  entrainment,  it  is  also  globally  limited  by  it.  The  mechanism  of 
this  apparent  double  role  is  explained  as  follows.  When  the  Reynolds  number 
is  high,  diffusion  between  the  reactants  occurs  across  a  highly  convoluted 
interfacial  layer  embedded  within  the  large  structures  at  a  much  faster  rate 
than  the  entrainment  of  the  reactants  into  the  large  eddies.  Thus,  changing 
the  diffusion  rate  by  varying  the  Reynolds  number  does  not  affect  the 
product  formation  since  the  process  is  entrainment  limited. 

The  phenomenon  of  mixing  asymmetry,  which  arises  due  to  the  finite 
velocity  difference  between  the  two  streams  [7],  is  shown  in  the  product 
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distribution  in  figure  14b,  and  is  manifested  in  the  dependence  of  the 
product  formation  of  the  reactants  ratios  in  figure  15c.  The  reactants 
ratio  is  defined  as  the  ratio  between  the  two  reactants  concentrations 
across  the  layer.  For  equal  concentrations  of  reactants  in  the  two  streams, 
the  average  product  concentration  in  the  lower,  slower  stream  is 
persistently  higher  them  in  the  upper,  faster  stream.  This  is  because  the 
large  structures  entrain  more  high-speed  them  low-speed  fluid.  The  bias 
towards  the  entrainment  of  the  high-speed  stream  results  in  a  faster 
chemical  reaction  in  the  lower  section  of  the  layer  than  in  the  upper 
section.  Increasing  the  concentration  of  either  reactants  enhances  the  rate 
of  product  formation  since  the  chemical  reaction  becomes  faster.  However, 
to  improve  the  overall  rate  of  burning,  it  is  more  effective  to  increase  the 
slow-stream  concentration  to  compensate  for  the  entrainment  asymmetry. 

VI. 2.  3D  reacting  shear  layers. 

We  have  shown  in  Section  III. 3  that  three-dimensional  secondary 
instabilities  follow  the  saturation  of  the  two-dimensional  instabilities. 
The  complicated  mixing  patterns,  which  arise  due  to  the  evolution  of  three- 
dimensional  instabilities  in  shear  layers,  are  revealed  by  the  results  of  a 
reacting  flow  simulation  between  a  stream  of  fuel  and  a  stream  of  oxidizer 
[30].  Figure  16  shows  product  concentration  and  reaction  rate  contours  on 
(a)  a  spanwise  x-z  section  which  passes  through  the  midsection  of  the  domain 
of  figure  7;  and  (b)  and  (c)  the  two  streamwise  y-z  sections  used  to  depict 
vorticity  in  figures  8  and  9.  An  incompressible,  temporal,  periodic  flow 
model,  and  a  reaction  rate  W  =  D  c c_,  are  used.  The  physical  parameters 
are:  D_  *  10,  P  *  1000  and  L  =1.  We  only  show  the  late  staqes  when  three 
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dimensional  effects  become  important. 


35 


These  figures  show  an  unexpected  inhomogeneity  of  mixing,  and 
combustion  inside  the  cores  in  the  streamwise  and  the  spanwise  directions, 
even  at  the  late  stages  of  instability  development.  As  mentioned  before, 
these  flows  are  entrainment  dominated.  The  mixing  between  the  two  streams 
is  initiated  along  convoluted  interfaces  which  are  produced  by  the  uneven 
stretch  within  the  large  structures.  The  stretch  is  a  consequence  of  the 
vorticity  instability.  Diffusion,  a  process  much  slower  than  convection  at 
high  Reynolds  number,  occurs  within  these  thin,  convoluted  zones.  The 
chemical  reaction  at  D  -  10  converts  mixed  reactants  into  products  almost 
immediately.  Thus,  products  are  formed  within  thin  convoluted  zones 
embedded  inside  the  large  structures. 

In  two-dimensional  simulations,  products  are  found  inside  the  large 
cores.  In  three-dimensional  simulations,  products  are  found  within  the 
large  cores  and  within  the  braids  due  the  formation  of  the  rods.  In  both 
cases,  zones  of  maximum  strain  in  the  spanwise  and  streamwise  planes  are 
depleted  of  products  by  mechanisms  which  will  be  discussed  next  section. 
Thus,  3D  simulation  confirm  our  earlier  observation  that  product 
concentration  is  maximum  where  vorticity  concentration  is  highest,  and  that 
at  zones  of  high  strain  and  small  vorticity,  product  concentration  is 
lowest.  This  correlation  holds  at  all  times  and  across  all  sections.  This 
correlation  has  been  observed  in  the  spanwise  plane  in  several  two- 
dimensional  computations  [40].  In  the  streamwise  planes,  the  correlation 
between  vorticity  and  product  concentration  is  supported  by  figures  16b  and 
c;  and  figures  8  and  9. 

The  dynamic  origin  of  this  correlation  is  revealed  by  inspecting  the 
reaction  rate  contours  shown  in  figure  17.  At  the  early  stages,  products 
form  at  the  center  of  the  eddies  where  mixing  occurs  due  to  the  entrainment 
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field.  Later,  reaction  is  maximum  at  the  outer  edges  of  the  structure. 
Thus,  at  this  high  value  of  Da,  product  formation  zones  move  from  the  center 
of  the  core  to  the  outer  edges  of  the  eddies.  However,  at  all  stages,  the 
entrainment  field  draws  the  products  towards  zones  of  maximum  vorticity 
which  exist  inside  the  cores  and  the  rods.  Thus,  although  products  may  form 
away  from  the  zone  of  maximum  vorticity,  they  are  constantly  brought  there 
by  the  convective  field. 

The  reaction  zone  does  not  exhibit  a  similar  correlation  with  the 
vorticity.  It  has,  thus,  been  suggested  that  vorticity  be  used  as  a  primary 
variable  in  turbulent  combustion  models  [41]. 

VI. 3.  THE  REACTING  JET. 

The  effect  of  exothermic  energy  on  the  flow  dynamics  and  the  structure 
of  the  reaction  zone  was  investigated  in  the  case  of  a  reacting  jet  of  fuel 
issuing  in  an  atmosphere  of  oxidizer  [27,42].  The  flow  is  two-dimensional 
and  planar  and  the  jet  forms  two  vorticity  layers  with  two  opposite  signs  on 
both  sides  of  its  centerline.  We  used  a  compressible  flow  model,  Eqs.  (21)- 
(26),  and  a  single  step  Arrhenius  reaction,  w  »  Af  c_  c  exp  (-T /T). 
Figure  17  depicts  the  product  distribution  for  Af  -  150  and  750.  Both  were 
obtained  at  Pg  *  1000,  Lg  -  1,  Tp/TR  *  4,  Q  -  3,  and  Ta/TR  -  10.  The 
initial  perturbation  amplitude  is  e  *  0.05  X,  where  X  is  the  wavelength  in 
the  streamwise  direction.  Figure  18  shows  the  reaction  rate  contours  for 
the  two  cases. 

VI. 3.1.  PRODUCT  DISTRIBUTION 

At  high  Damkohler  number,  product  concentration  is  higher  due  to  the 
faster  chemical  reaction,  and  the  eddies  are  larger  due  to  the  extra  heat 
release,  as  shown  in  figure  16.  At  both  values  of  D  ,  product  concentration 
within  the  shear  layers  changes  from  a  uniform  distribution  in  the 
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streamwise  direction  to  a  highly-concentrated  distribution  within  the  large 
eddies  as  the  flow  evolves.  Product  concentration  within  the  braids, 
however,  decreases  continuously.  The  thinning  of  the  braids,  due  to  the 
formation  of  a  strong  strain  field  as  the  eddies  roll-up,  reduces  the  zones 
of  overlap  between  the  jet  and  oxidizer  fluids,  of  the  zone  of  chemical 
activity.  Moreover,  the  large  relative  velocity  between  the  jet  and 
oxidizer  streams  within  the  braids  also  reduces  the  time  available  for 
mixing  within  this  region.  Finally,  the  entrainment  currents  convects 
products  forming  within  the  braids  into  the  cores  leading  to  continuous 
cooling  of  the  region  between  neighboring  eddies.  This  the  physical 
mechanism  of  reaction  extinction  within  the  braids.  Clearly,  it  is 
convection-driven. 

The  similarity  between  product  distribution  and  vorticity,  observed  in 
Section  VT.l,  persists  in  the  case  of  a  reacting  jet,  both  at  low  and  high 
Damkohler  number.  Both  product  concentration  and  vorticity  exhibit  high 
values  within  the  eddies  and  fall  continuously  along  the  braids.  Thus, 
vorticity  still  controls  product  distribution  within  the  eddy.  At  high 
Damkohler  number,  vorticity  produces  a  swirling  field  that  enlarges  the 
surface  of  contact  between  the  two  reactants,  diffusion  mixes  them  across 
the  stretched  surface,  and  products  form.  The  swirling  convective  field 
then  entrains  these  products  into  the  eddy  core.  At  low  damkohler  number, 
reactants  are  entrained  into  the  cores  and  then  react.  The  difference 
between  the  two  mechanisms  is  illustrated  by  examining  the  reaction  rate 
contours  for  both  cases. 

VI. 3.1.  STRUCTURE  OF  THE  REACTION  ZONE. 

At  low  Damkohler  number,  reaction  rate  contours  presented  in  figure  17 
show  that  the  zone  of  chemical  activity,  initially  uniform  in  the  streamwise 
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direction,  is  enlarged  during  the  entrainment  phase  of  the  eddy.  The  mixing 
of  jet  and  ambient  fluids,  induced  by  the  rotational  field  within  the  eddy 
core,  extends  the  area  of  combustible  mixture.  Moreover,  the  reaction  rate 
inside  the  eddies  increases  with  time  triggered  by  the  acceleration  of  the 
chemical  reaction  as  the  temperature  rises  due  to  the  formation  of  more 
products.  Large  eddies,  thus,  act  as  exothermic  centers  that  support 
combustion.  At  this  low  Damkohler  number,  an  eddy  acts  as  "mixture 
preparation  zone"  before  combustion  proceeds.  This  is  consistent  with  the 
fact  that  the  rotational  flow  within  the  eddy  occurs  at  a  smaller  time  scale 
than  the  chemical  reaction. 

At  high  Damkohler  number,  the  reaction  zone  differs  substantially  from 
that  at  low  Damkohler  number.  At  early  stages,  the  reaction  rate  is  maximum 
inside  the  eddy  where  combustion  occurs  over  a  distributed  zone.  In  this 
case,  both  reactants  are  drawn  into  the  mixing  zone  within  the  eddy  core 
before  they  combust.  At  late  stages,  combustion  occurs  on  the  outer  edges 
of  the  eddies  within  a  thin  reaction  zone.  Now,  reactants  coexist  only 
around  the  outer  edges  of  each  eddy  where  they  react  and  form  products  which 
are  then  entrained  into  the  eddy  core.  Note  that  with  finite-rate  kinetics, 
and  with  low-temperature  reactants  at  the  early  stages,  there  is  an  ignition 
delay  which  keeps  the  kinetic  rate  lower  than  the  mixing  rate.  As  the 
temperature  rises,  kinetic  rates  exceed  mixing  rates  and  the  reaction  region 
becomes  a  thin  zone  between  jet  and  ambient  fluids. 

At  low  and  high  D  ,  the  reaction  rate  within  the  braids  drops  sharply 
as  the  eddies  form  due  to  the  entrainment  of  products  into  the  eddy  cores 
(Note  that  since  the  Lewis  number  is  unity,  temperature  and  piodu'-t 
concentration  are  similar).  Product  formation  is  inhibited  within  the 
braids  following  their  cooling  by  convection. 
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VI. 3. 3.  EFFECT  OF  HEAT  RELEASE. 

In  a  reacting  flow  with  finite-rate  kinetics,  we  found  that  heat 
release  reduces  that  rate  of  growth  of  the  instability  [27]  when  the  initial 
perturbation  is  small,  e  -  0.01  X.  The  suppression  of  the  instability 
becomes  stronger  as  the  Damkohler  number  is  increased.  Similar  suppression 
of  instability  is  observed  at  infinite-rate  kinetics  at  a  range  of  initial 
perturbation  e  *  0.01  -  0.05  X  [42].  In  both  cases,  eddy  growth  in  the 
direction  normal  to  the  jet  axis  is  reduced  as  the  Damkohler  number  is 
increased.  The  reason  is  as  follows:  At  infinite-rate  kinetics,  exothermic 
energy  and  volumetric  expansion  start  early,  within  the  linear  stages  of 
instability  growth,  forming  a  sublayer  of  low-density  products  within  the 
vorticity  layer.  Within  these  stages,  it  has  been  shown  that  the  presence 
of  low-density  sublayer  within  the  vorticity  layer  leads  to  substantial 
reduction  of  the  instability  growth  rate  [43].  The  delay  in  roll-up  results 
in  the  formation  of  a  weak,  less  coherent  eddy.  Even  at  finite-rate 
kinetics  with  small  initial  perturbation,  an  appreciable  amount  of  products 
forms  within  the  vorticity  layer  before  the  instability  amplitude  is  large 
enough  to  initiate  roll-up.  In  this  case,  the  early  formation  of  a  low- 
density  sublayer  of  products,  embedded  within  the  vorticity  layer,  precedes 
the  rollup  and  leads  to  an  overall  reduction  in  the  instability  growth  rate. 

At  finite-rate  kinetic  with  large  initial  amplitude,  e  -  0.01  X, 
instability  suppression  is  negligible,  as  seen  in  Figure  17.  Thus,  one  can 
overcome  the  stabilizing  effect  of  volumetric  expansion  on  shear  layer 
growth  by  forcing  the  jet  at  large  amplitudes.  Initiating  roll-up,  before 
measurable  volumetric  expansion  has  occurred,  causes  eddy  formation  and  the 
onset  of  mixing  enhancement  during  the  ignition  delay  time.  Very  small 
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initial  perturbations  with  high  exothermic  energy,  on  the  other  hand,  result 
in  instability  suppression. 

VI. 4.  THE  PREMIXED  SHEAR  LAYER 

Computations  of  a  shear  layer  growing  between  two  streams  of  premixed 
reactants  and  products  were  presented  in  [7,35].  In  this  case,  the  reaction 
rate,  based  on  one-step  reaction  kinetics  R  ->  P,  is  expressed  as  W  ■  A^  cR 
exp(-T/T  ).  Results  for  a  range  of  Damkohler  number  and  initial  amplitude 
of  perturbations  are  presented.  A  sample  of  the  results  is  shown  in  figure 
19  in  terms  of  the  rate  of  reaction,  the  product  concentration,  and  the 
vorticity  for  the  case  of  Ta/TR  -  10,  temperature  ratio  Tp/TR-  5,  Q  -  4,  A^- 
1,  Pg  «  1000  and  L  •  1.  Figure  20  shows  the  total  mass  of  products  formed 
within  the  layer  for  the  same  physical  parameters  but  for  several  values  of 
the  Damkohler  number.  The  length  of  the  line  of  maximum  reaction  rate,  and 
the  total  mass  of  product  for  a  representative  case  is  shown  in  figure  21. 
A  careful  inspection  of  these  results  reveal  several  interesting 
observations. 

Figures  19  and  20  indicate  that  one  can  divide  the  combustion  in  a 
premixed  shear  layer  into  four  phases:  (1)  a  laminar  flame;  (2)  a  strained 
laminar  flame;  (3)  a  vortex-driven  combustion,  (4)  a  free-propagating  flame. 
In  the  second  phase,  the  reaction  zone  in  figure  19  is  thinner  than  that  of 
a  laminar  flame  due  to  the  strain,  and  the  total  amount  of  products  formed, 
Mp  *  J  p  cp  dy  dx,  is  less  for  the  shear  layer  them  for  the  laminar  flame, 
Mp  «  X  where  X  is  the  perturbation  wavelength.  In  the  third  phase,  the 
entrainment  associated  with  the  formation  of  a  coherent  vortex  leads  to  the 
swelling  of  the  core  and  the  establishment  of  a  reaction  zone  inside  t he 
eddy  core.  Figure  20  shows  that  the  formation  of  a  vortex  core  enhances  the 
rate  of  burning.  As  the  flame  leaves  the  burning  vortex,  it  returns  to  the 
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state  of  a  laminar  flame.  Since  these  phases  result  from  the  interaction 
between  the  flow  and  the  chemical  reaction,  their  starting  time  and  duration 
are  strongly  dependent  on  the  Damkohler  number. 

Figure  21  shows  that  the  wrinkled  flame  model,  which  states  that  the 
rate  of  product  formation  -  Su  where  Su  is  the  laminar  burning 
velocity  and  is  the  total  flame  length,  earn  be  used  to  approximate  the 
burning  rate  during  the  initial  growth  phase  of  the  eddy  provided  that  Lf  is 
measured  along  the  line  of  maximum  reaction  rate.  However,  during  the  later 
stages,  the  value  of  Su,  as  defined  above,  is  found  to  decrease  below  the 
value  of  the  unstrained  flame.  During  the  later  stages,  after  most  of  the 
eddy  core  has  burnt,  it  is  difficult  to  define  a  flame  front  due  to  the 
convolution  of  the  streamlines  and  further  studies  of  the  detail  structure 
of  the  reaction  zone  is  needed. 

Baroclinic  vorticity  generation,  associated  with  the  interaction 
between  the  hydrodynamic  pressure  gradient  and  density  field,  contributes  to 
the  dynamics  of  the  shear  layer  in  the  same  way  as  in  the  nonreacting, 
density-stratified  shear  layer  described  earlier:  the  large  eddies  move  in 
the  direction  of  the  cold  reactants;  entrainment  asymmetry  biases  the 
composition  of  the  large  eddies  towards  the  hot  product  stream;  and  local 
spotiness  appears  within  the  large  eddy.  One  difference  is  clear  by 
comparing  figures  12  and  19:  heat  release  slows  down  the  streamwise  motion 
of  the  large  eddies.  Due  to  the  volumetric  dilatation  associated  with  heat 
release,  the  local  concentration  of  vorticity  decreases  and  the  induced 
field  on  the  vortex  weaken.  Thus,  heat  release  weaken  the  instability 
through  the  mechanism  of  volumetric  expansion  and  not  vorticity  generation. 

Heat  release  at  high  Damkohler  number  reduces  the  growth  rate  of  the 
instability  at  the  initial  stages  and  inhibits  the  rollup  process  if  the 
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initial  perturbation  is  small,  e  =  0.01  X,  similar  to  the  reacting  jet  case. 
The  mechanism  of  instability  suppression  in  this  case  is  different  from  that 
found  in  the  reacting  shear  layer.  In  a  premixed  shear  layer,  baroclinic 
vorticity  by  itself  enhances  the  instability  growth  rates,  as  shown  before. 
Volumetric  expansion  is,  thus,  necessary  for  instability  suppression  in  the 
premixed  shear  layer  flow. 

The  effect  of  the  Damkohler  number  on  the  structure  of  the  reaction 
zone  is  similar  to  the  reacting  jet;  at  low  Damkohler  number,  reaction 
occurs  within  the  eddy  core  and  at  high  Damkohler  number,  it  occurs 
primarily  on  the  outer  edges  of  the  eddy.  As  in  the  reacting  shear  layer, 
product  concentration  is  always  higher  at  the  center  of  the  eddy  (premixed 
shear  layer  are  however  more  subtle  since  product  forming  during  the 
development  of  the  shear  layer  and  products  existing  at  the  initial  state 
ate  not  easily  distinguishable).  As  shown  in  figure  20,  the  effect  of  the 
shear  layer  on  the  burning  rate  is  stronger  at  low  Damkohler  number. 

VII.  EXTENSIONS 

Work  reviewed  so  far  focuses  on  combustion  in  free  shear  flows.  Brief 
summary  of  the  numerical  methodology  and  its  application  to  nonreacting  and 
reacting  shear  layer  has  been  presented.  Analysis  of  the  computational 
results  aimed  at  testing  the  accuracy  of  the  numerical  schemes  and 
describing  some  of  the  general  properties  of  these  flows.  In  the  reacting 
flow,  mechanisms  of  shear  flow-combustion  interactions  were  analyzed  in 
light  of  the  numerical  solutions. 

The  application  of  vortex  methods  to  internal,  wall  bounded  flows  in 
which  the  growth  of  boundary  layers  along  solid  walls  and  their  separation 
at  sharp  edges  play  a  dominant  role  in  the  dynamics  of  the  flow,  has  been 
based  largely  on  the  random  vortex  method.  In  this  method,  the  effect  of 
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molecular  viscosity  is  taken  into  account  by  adding  a  Gaussian  random 
component  to  the  convective  motion  of  the  vortex  elements  [44,45]. 
Extensive  work  on  the  validation  of  the  method  [46,47]  shows  that  solutions 
for  steady,  low  Reynolds  number  flows,  and  unsteady,  high  Reynolds  number 
flows  converge  to  appropriate  limits  as  numerical  parameters  are  refined. 
Low-Reynolds-number  results  were  in  agreement  with  experimental  measurements 
on  velocity  distributions  within  the  flow.  At  high  Reynolds  number,  results 
were  shown  to  converge  to  oscillatory  flows  that  can  be  characterized  by 
time-dependent  clusters  of  large  scale  vortices.  The  dependence  of  the 
shedding  frequency  and  growth  rate  of  these  structures  on  the  geometry  were 
investigated  in  [48].  The  interactions  between  jet  flow,  recirculating  flow 
and  annulus  flow,  encountered  in  bluff-body  flame  burner,  has  been  revealed 
in  a  vortex  simulation  [49]. 

The  random  vortex  method  was  also  applied  to  study  reacting 
recirculating  flows  of  premixed  gases  at  high  Reynolds  numbers  [50,51] 
utilizing  the  thin  flame  approximation  [52]  to  model  the  combustion  process. 
Results  were  used  to  study  vorticity-flame-pressure  interactions  in  a 
semiconfined,  recirculation-stabilized  premixed  flame.  Analysis  of  these 
results  show  how  the  volumetric  expansion  associated  with  burning  can  reduce 
the  amplification  of  flow  oscillation  in  recirculating  flows  when  the 
pressure  is  kept  constant,  and  how  these  oscillation  cam  be  amplified 
leading  to  large  flame  oscillation  if  the  pressure  field  within  the  system 
is  coupled  with  the  flow  processes  [51,53].  These  results  address  the 
problem  of  turbulent  flame  stabilization  in  a  premixed  stream  and  the 
associated  combustor  instabilities  observed  in  such  systems. 


VIII.  CLOSURE 
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Numerical  simulation,  using  accurate  schemes  to  integrate  the  unsteady 
equations  governing  reacting  flow,  can  be  applied  to  investigate  important 
mechanisms  of  shear  flow-combustion  interaction  in  different  systems  and 
within  a  wide  range  of  physical  parameters.  Mechanisms  of  flow  combustion 
interaction  are  different  in  nonpremixed  and  premixed  shear  layers. 
Entrainment,  associated  with  the  formation  of  structures  that  develop  due  to 
natural  shear  flow  instability,  is  the  primary  mechanism  by  which  the  flow 
enhances  the  rate  of  burning  in  both  cases.  Burning  enhancement  is 
substantial ,  especially  at  low  Damkohler  number  where  entrainment  changes 
the  reaction  region  from  a  thin  front  into  a  distributed  zone.  Strain  may 
have  some  effect  on  the  burning  mechanisms  especially  at  low  Damkohler 
number  where  local  extinction  has  been  observed  at  regions  of  high  strain 
rate  in  the  nonpremixed  shear  layer.  In  the  premixed  shear  layer,  the 
reduction  in  the  rate  of  burning  due  to  strain  is  greater  at  high-Damkohler 
number. 

Heat  release  establishes  zones  of  density  gradient  within  the  vorticity 
layer.  Baroclinic  vorticity,  generated  from  the  interaction  between  these 
density  gradients  and  material  acceleration,  reduces  the  growth  rate  of  the 
instability  in  the  nonpremixed  case  while  enhancing  it  in  the  premixed  case. 
The  later  occurs  even  in  the  nonreacting  case.  Volumetric  expansion, 
however,  suppresses  the  instability  in  both  cases.  The  effect  is  more 
pronounced  at  high-Damkohler  number.  Volumetric  dilatation  reduces  the 
vorticity  and  weaken  the  large  structures.  Forcing  at  high  amplitudes  in 
the  initial  stages  can  be  used  to  overcome  this  instability  suppression. 

Work  is  under  way  to  extend  the  numerical  methodology  in  different 
directions,  e.g.,  (1)  the  formulation  of  a  transport  element  method  for 
confined  flows  in  which  the  interaction  between  the  interior  flow  and  the 
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wall  thermal  boundary  layer  is  important;  (2)  the  application  of  the 

transport  element  method  to  flows  in  which  pressure-density  interaction 

plays  3  significant  dynamic  rol.-  [54];  (3)  the  implementation  oi  fast 

2 

solvers  to  reduce  the  computational  effort  from  0(N  )  to  0(N)  [55];  (4)  the 
implementation  of  a  methodology  to  allow  the  application  of  the  scheme  to 
extended,  multistep  chemical  kinetic  models;  and  (5)  the  formulation  of 
Lagrangian  schemes  for  multiphase  flow  such  as  fuel  droplets  in  gas  streams. 
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FIGURE  CAPTIONS 

Figure  1.  The  rollup  of  an  incompressible,  uniform  density,  nonreacting, 
temporally-growing  shear  layer  in  which  the  top  and  bottom  streams  move  at 
the  same  speed  but  in  opposite  directions  and  the  initial  vorticity 
distribution  is  guussian.  The  initial  amplitude  of  the  sinewave 
perturbation  is  1%  of  its  wavelength,  which  is  taken  as  the  most  unstable 
wavelength  for  this  flow.  The  results  are  shown  in  terms  of  (a)  the  vortex 
elements  and  their  velocity  vectors;  and  (b)  the  isoscalar  lines  (s  -  0  in 
the  top  stream  and  s  -  1  in  the  bottom  stream  with  an  error  function 
distribution  in  between). 

Figure  2.  The  development  of  an  incompressible,  nonreacting,  spatially 
developing  shear  layer  between  a  fast,  top  stream  and  a  slow,  bottom  stream. 
The  velocity  ratio  across  the  layer  is  2  and  the  vorticity  distribution  at 
the  left-hand  side  of  the  domain,  where  the  splitter  plate  ends,  is 
gaussian.  The  figures  show  all  the  vortex  elements  used  in  the  computation 
and  their  velocity  vector  (measured  with  respect  to  the  mean  velocity  of  the 
flow).  The  confining  walls  are  slip  boundaries.  Three  cases  are  shown:  (a) 
the  layer  is  unforced.  (b)  the  layer  is  forced  at  the  most  unstable 
wavelength  and  the  amplitude  of  the  perturbation  is  1%  of  the  wavelength, 
(c)  The  layer  is  forced  at  the  most  unstable  wavelength  and  its  first 
subharmonic,  and  the  two  amplitudes  are  equal  to  1%  of  the  fundamental 
wavelength. 

Figure  3.  (a)  Time-average  streamwise  velocity  profiles  for  the  case  shown 
in  figure  2a,  computed  at  sections  x-3,  3.5,  4,  4.5  and  5.  x  is  measured 
from  the  left-hand  side  of  the  domain,  x  -  xO,  and  is  normalized  with 
respect  to  the  channel  height,  while  y0-0.  Uh  and  Ul  are  the  high  and  low 
speed  stream  velocities.  Computational  results  are  shown  in  solid  lines  and 
experimental  measurements  of  Masutani  and  Bowman  [29]  are  shown  in  open 
symbols. 

(b)  Time-average  streamwise  velocity  fluctuations  for  the  same  case  as  in 
(a),  computed  at  the  same  sections  and  plotted  against  experimental 
measurements  from  the  same  reference. 

Figure  4.  The  deformation  of  a  thin  vortex  ring,  represented  computationally 
by  a  single  vortex  filament,  when  excited  at  an  unstable  wavenumber  (number 
of  waves  around  the  ring)  of  six.  The  radius  of  the  ring  core  is  0.25  of 
the  radius  of  the  ring.  The  plots  are  obtained  by  projecting  the  ring  on 
two  planes  normal  (on  the  left-hand  side)  and  parallel  (on  the  right-hand 
side)  to  the  direction  of  propagation  of  the  ring  under  its  own  self-induced 
velocity. 

Figure  5.  The  wavenumber  of  the  most  unstable  modeT  n*,  of  a  vortex  ring 
plotted  against  its  normalized  self-induced  velocity,  V  =  4nRV/T  where  R  and 
T  are  the  radius  and  circulation  of  the  ring,  respectively.  The  figure 
shows  a  comparison  between  the  experimental  results,  x;  the  analytical 
results  of  the  long  wavelength  instability,  o;  the  numerical  results  of  the 
long  wavelength  instability,  A;  the  analytical  results  of  the  short 
wavelength  instability,  for  constant  vorticity  distribution  within  the 
core  and  +  for  quadratic  vorticity  distribution;  and  the  computed  results  of 
the  short  wavelength  instability,  <>  for  coarse  numerical  discretization  and 
*  for  fine  discretization. 


Figure  6.  Perspective  views,  taken  form  the  point  of  view  of  an  observer 
standing  ahead  of  the  ring  and  looking  at  an  angle  60  with  respect  to  the 
direction  of  propagation,  of  a  vortex  ring  whose  core  radius  is  0.275  of  the 
ring  radius.  The  ring  is  initially  excited  using  12  waves  around.  All  the 
filaments  used  to  discretize  the  ring  are  shown.  The  ring  is  propagating 
upwards. 

Figure  7.  Three-dimensional  perspective  views  of  the  isoscalar  surface  s  - 
0,  initially  coinciding  with  an  x-y  plane  located  in  the  middle  of  the 
domain  shown  in  the  figure.  The  shear  layer  is  periodic  in  two  directions. 
At  time  t  *  0,  the  layer  is  perturbed  in  the  streamwise  direction,  x,  at  the 
most  unstable  wavelength,  and  in  the  spanwise  direction,  y,  at  one-half  of 
this  wavelength.  The  amplitudes  of  the  perturbations  are  equal  to  2%  of  the 
streamwise  wavelength.  The  initial  vorticity  follows  a  gaussian 
distribution  in  the  spanwise  direction,  z,  as  in  the  two-dimensional  case. 

Figure  8.  (a)  Streamwise  vorticity;  and  (b)  scalar  distribution  in  the 
three-dimensional  simulation  of  the  periodic  shear  layer  shown  in  figure  7. 
Contours  are  shown  in  a  y-z  plane  which  cuts  through  the  core  of  the 
spanwise  eddy  at  the  middle  of  the  domain  of  figure  7.  Positive  vorticity 
is  indicated  by  broken  lines  while  negative  vorticity  is  depicted  by 
continuous  lines.  For  the  scalar,  s  >  0  is  shown  in  continuous  lines  and  s 
<  0  is  shown  in  broken  lines,  -0.5  <  s  <  0.5. 

Figure  9.  (a)  Streamwise  vorticity;  and  (b)  scalar  distribution  in  the 
three-dimensional  simulation  of  a  periodic  shear  layer  shown  in  figure  7. 
Contours  are  shown  in  the  y-z  plane  which  cuts  through  the  braids  at  a 
distance  15%  of  the  wavelength  into  the  domain  from  the  left-hand  side.  See 
figure  8  for  convention. 

Figure  10.  (a)  Time-average  scalar  concentration  for  the  flow  in  figure  2a, 
computed  at  x  -  5  for  a  range  of  Peclet  number,  P  .  At  x  ■  0,  the 
concentration  is  s  -  0  in  the  upper  stream  and  s  ■  1  in  the  lower  stream. 
Solid  lines  show  the  results  of  the  computations  and  open  symbols  depict  the 
experimental  measurements  of  Masutani  and  Bowman  [29]. 

(b)  Time-average  concentration  fluctuations  for  the  same  flow  as  in  (a). 

Figure  11.  The  eddy  size  in  the  two,  v;  and  three  ,  x;  dimensional 
computations  of  figures  1  and  7,  respectively.  The  eddy  size  is  defined  by 
the  cross-stream  distance  between  the  contours  s  -  0.03  and  s  -  0.97. 

Figure  12.  The  rollup  of  an  incompressible,  nonuniform  density,  nonreacting, 
temporally-growing  shear  layer  in  which  the  top,  heavy  stream  moves  to  the 
right  and  the  bottom,  light  stream  moves  to  the  left  at  the  same  speed.  The 
initial  conditions  are  the  same  as  in  figure  1,  except  for  the  density 
distribution  across  the  layer  which  is  taken  as  an  error  function.  1  <  p  < 
5.  The  density  ratio  across  the  layer  is  five.  The  results  are  shown  in 
terms  of:  (a)  the  isoscalar  lines  (s  =  0  and  in  the  top  stream  and  s  =  1  in 
the  bottom  stream);  and  (b)  the  vorticity  contours.  Broken  lines  indioat0 
positive  vorticity  and  solid  lines  indicate  negative  vorticity. 

Figure  13.  Three-dimensional  perspective  views  of  the  isoscalar  surface  s  = 
0,  initially  coinciding  with  an  x-y  plane  in  the  middle  of  the  domain.  The 
flow  is  incompressible  but  with  variable  density.  Initially,  the  streamwise 
velocity  and  density  vary  only  in  the  cross-stream,  z-direction  according  to 
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similar  error  functions  with  different  boundary  conditions:  -1  <  u  <  1  and  1 
<  p  <  2.  The  shear  layer  is  periodic  in  the  x  and  y  directions.  All  other 
conditions  are  similar  to  the  case  presented  in  figure  7. 

Figure  14.  Results  of  a  computation  of  a  spatially-developing,  reacting 
shear  layer  using  an  incompressible  flow  model  and  Arrhenius  chemical 
kinetics.  Figure  2a  shows  the  vortex  elements  are  their  velocity  vectors 
for  this  case  (since  the  flow  is  incompressible,  the  reaction  does  not 
affect  the  dynamics).  The  figure  shows  contours  of  (a)  the  top  stream 
reactant's  concentration;  (b)  product  concentration;  (c)  vorticity;  and  (d) 
temperature. 

Figure  15.  Results  of  a  spatially-developing  reacting  shear  layer  using  an 
incompressible  flow  model  and  temperature-independent  kinetics.  The  dynamic 
field  is  the  same  as  that  shown  in  figure  2a.  The  figure  shows  (a)  the 
product  integral,  /  e  dy,  across  the  layer  for  a  reacting  perturbed  shear 
layer,  similar  that  shown  figures  2  and  14,  and  across  an  unperturbed  shear 
layer  where  rollup  is  suppressed;  (b)  the  product  integral  as  a  function  of 
the  Damkohler  number;  (c)  The  product  integral  as  a  function  of  the  Reynolds 
number,  the  lines  are  computed  for  R  -  5000  and  10000  while  the  points  1,  2 
,  and  3  are  computed  at  R  *  500,  2500  and  50000,  respectively,  D  is  a 
Damkohler  number  based  on  the  local  value  of  x;  and  (d)  the  total  amSfint  of 
product  in  the  field,  f  f  c_  dy  dx  as  a  function  of  the  reactants  ratio 
across  the  layer. 

Figure  16.  Results  of  a  three-dimensional  simulation  of  doubly  periodic 

reacting  shear  layer  between  two  reacting  stream.  The  dynamics  of  this  flow 

is  the  same  as  that  in  figure  7  since  the  flow  is  incompressible.  The 

figure  shows:  (a)  the  product  concentration  (left)  and  reaction  rate  (right) 

on  an  x-z  spanwise  plane  located  at  the  middle  of  the  domain  in  figure  7. 

(b)  The  product  concentration  (left)  and  reaction  rate  (right)  across  the 

section  described  in  figure  8.  (c)  The  product  concentration  (left)  and 

reaction  rate  (right)  across  the  section  described  in  figure  9.  Results  are 

shown  in  gray  scale  in  which  the  maximum  value  is  indicated  by  black  and  the 

minimum  value  is  shown  in  white.  The  gray  scale  for  the  product 

concentration  is  fixed,  0  <  c_  <  1,  while  the  gray  scale  for  the  reaction 

rate  is  floating,  W  . _  <  W  <  w 

min  -  -  max 

Figure  17.  Product  distribution  within  the  initial  stages  in  a  reaction  jet 
using  a  compressible  flow  model  with  Arrhenius  chemical  kinetics.  The  model 
is  periodic  in  the  streamwise  direction.  Results  are  shown  at  three  time 
steps  for:  (a)  A^  »  150;  and  (b)  Af  -  750. 

Figure  18.  The  distribution  of  the  reaction  rate  within  one  side  of  the 
centerline  of  the  reacting  jet  of  figure  16,  shown  using  a  gray  scale  for 
(a)  Af  -  150  and  (b)  A,  -  750.  Results  are  shown  in  gray  scale  in  which  the 
maximum  value  is  indicated  by  black  and  the  minimum  value  is  shown  in  white. 

Figure  19.  The  roll-up  of  a  compressible,  temporally  growing  reacting  shear 
layer  between  cold  premixed  reactants  in  the  top  stream  and  hot  products  in 
the  bottom  stream.  Results  are  shown  in  terms  of:  (a)  the  reaction  rate 
(shown  in  gray  scale);  and  (b)  the  product  concentration  contours.  The 
temperature  ratio  across  the  layer  is  five. 
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Figure  20.  The  total  mass  of  products  formed  within  the  premixed  reacting 
shear  layer  shown  in  figure  19,  compared  with  the  same  quantity  for  a 
laminar  flame  propagating  through  the  same  mixture.  Results  are  shown  for 
three  different  values  of  Af.  Straight  lines  show  the  total  mass  of 
products  formed  in  the  corresponding  laminar  flame. 

Figure  21.  The  total  mass  of  products  and  the  length  of  the  line  of  maximum 
reaction  rate  in  a  premixed  reacting  shear  layer  similar  to  the  one  shown  in 
figure  19,  and  in  the  corresponding  laminar  flame.  The  temperature  ratio 
across  the  layer  is  three,  ■  200  and  T  =  5. 
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